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An  analytical  study  of  the  thermodynamic  processes  occurring 
in  an  Otto  cycle  engine  has  been  conducted  with  a  view  to  obtaining  a 
better  understanding  of  the  combustion  process  and  a  realistic 
prediction  of  the  emission  levels.   Much  of  the  work  previously 
reported  in  this  area  has  been  experimental,  the  reason  being  the 
extreme  complexity  of  any  realistic  solution  to  the  governing  equa- 
tions.  The  availability  of  high-speed  computing  machines  provides  the 
necessary  tool  for  such  an  approach.   With  such  machines,  theoretical 
calculations  have  assumed  a  new  importance  and  in  many  cases  are  being 
used  to  obtain  estimates  which  would  require  extensive  experimental 
work.   Theoretical  analysis,  by  its  very  nature,  involves  the  use  of 
simplifying  assumptions  to  reduce  any  real  problem  to  a  stage  where 
it  is  solvable.   Any  unrealistic  simplifying  assumption  will  reduce 
the  solution  to  one  of  academic  interest  only  and  the  literature  on 
any  subject  is  replete  with  such  examples.   In  this  study,  an  analysis 
with  a  minimum  of  sim.plifying  assumptions  has  been  carried  out  and  the 
results  obtained  are  realistic  and  usable. 
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The  study  presented  here  has  been  divided  into  two  parts. 
Firstly,  a  study  of  the  combustion  process  has  been  carried  out  by 
numerically  solving  the  set  of  partial  differential  equations  govern- 
ing one-dimensional  unsteady  laminar  flame  propagation  in  a  homogeneous 
fuel-air  mixture.   Secondly,  a  mathematical  model  for  the  Otto  cycle 
has  been  developed  to  predict  the  thermodynamic  properties  and  the 
species  concentrations  throughout  the  cycle  based  on  a  non-equilibrium 
analysis  and  finite  rate  combustion. 

For  the  first  part,  an  explicit  finite  difference  method  has 
been  used  to  numerically  solve  the  governing  equations  and  thereby 
obtain  a  detailed  description  of  the  flame  structure  in  terms  of  the 
distribution  of  the  state  properties  and  species  concentrations  through 
the  flame  zone.   Host  of  the  previous  analyses  have  ignored  the  varia- 
tions in  the  thermodynamic  and  transport  properties  with  temperature, 
pressure  and  chemical  composition  and  have  also  treated  the  chemical 
kinetics  as  being  governed  by  a  single  overall  reaction.   These 
particular  simplifying  assumptions  (particularly  the  last  one  which 
is  rather  unrealistic  for  any  system  of  practical  interest)  have  been 
removed  in  the  computational  scheme  presented  here.   The  scheme  is 
general  and  considers  variable  transport  and  thermodynamic  properties 
and  a  realistic  reaction  mechanism  in  terms  of  a  set  of  elementary 
reaction  steps.   The  validity  of  the  procedure  has  been  demonstrated 
by  carrying  out  calculations  for  the  methane-air  flame  and  obtaining 
good  agreement  with  the  limited  experimental  data  available. 

The  mathematical  model  of  the  Otto  cycle  gives  a  detailed 
description  of  the  cycle  in  terms  of  the  thermodynamic  properties 
and  species  concentrations  at  any  stage  in  the  cycle,  including  the 
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exhaust.   The  combustion  process  has  been  considered  to  proceed  at 
a  finite  rate  instead  of  the  commonly  made  assumption  that  this 
process  is  instantaneous.   Variations  in  the  thermodynamic  properties 
with  temperature  and  chemical  composition  have  been  allowed  instead 
of  assigning  them  constant  values.   An  accurate  prediction  of  the 
concentrations  of  the  various  species  at  exhaust  was  the  main 
objective  in  the  development  of  the  model. 

The  model  has  been  used  to  study  the  effect  of  the  intake 
of  an  additive  mixture  of  water  and  alcohols  (methanol  and  ethanol) 
in  different  proportions  and  amounts  on  the  emission  levels.   The 
model  is  quite  flexible  and  permits  the  study  of  the  influence  of 
other  factors  such  as  exhaust  gas  recirculation,  air-fuel  ratio, 
compression  ratio,  speed,  ignition  timing,  etc.,  on  the  performance 
and  emission  levels. 
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CHAPTER  I 
INTRODUCTION 

The  role  of  the  internal  combustion  engine  as  a  major 
contributor  of  the  hydrocarbon  gases,  carbon  monoxide,  and  oxides 
of  nitrogen  that  contaminate  the  atmosphere  has  been  recognized  for 
some  years.   The  formation  of  these  pollutants  is  closely  related 
to  the  coirbustion  process  and  a  realistic  modeling  of  this  process 
would  help  the  study  of  the  methods  to  reduce  the  emassions.   A 
knowledge  of  the  transient  mechanism  of  high  temperature  oxidation 
of  hydrocarbon  fuels  would  lead  to  a  better  understanding  of  the 
combustion  process  and  may  ultimatel}'  help  in  the  design  of  engines 
which  will  emit  lovjer  concentrations  of  pollutants.   I'Jliile  the  amount 
of  unbumt  or  partially  burnt  products  in  the  exhaust  is  generally 
insignificant  from  an  energy  point  of  view,  it  is  significant,  and 
very  much  so,  when  considered  from  an  environmentalist's  point  of 
view.   Thus  these  studies  leading  to  new  designs  might  or  might  not 
improve  the  overall  efficiency  of  the  internal  combustion  engine  (in 
fact  a  reduction  in  efficiency  might  be  achieved),  but  they  might 
help  to  reduce  the  emissions. 

The  high  temperature  oxidation  mechanism  of  some  simple  fuels 
like  hydrogen  and  carbon  monoxide  have  been  widely  studied  and  are 
fairly  well  understood,  but  the  state  of  development  for  the  oxidation 
of  hydrocarbon  fuels  leaves  much  to  be  desired.   Whatever  little  work 


has  been  reported  is  for  lower  hydrocarbons  (like  methane,  etc.)  and 
kinetic  data  and  reaction  mechanism  postulations  for  some  practical 
fuels  (like  iso-octane)  are  virtually  non-existent. 

One  of  the  methods  available  for  the  quantitative  determination 
of  chemical  kinetic  constants  is  a  detailed  analysis  of  the  flame 
structure.   A  comparison  of  the  flame  structure  predicted  theoretically 
on  the  basis  of  a  postulated  reaction  mechanism,  with  the  experimentally 
measured  data,  provides  a  convenient  \<iay   of  checking  the  validity  of 
the  reaction  mechanism.   The  influence  of  varying  the  rate  constant 
parameters  on  the  theoretical  predictions  permits  the  determination  of 
influence  coefficients  for  the  different  input  parameters.   A  parametric 
study  of  this  type  is  useful  in  arriving  at  some  sort  of  a  conclusion 
without  much  trial  and  error. 

To  put  this  principle  into  practice,  there  remains  the  need 
for  a  general  theoretical  computational  procedure,  and  it  is  the  aim 
of  this  study  to  proceed  in  this  direction  by  presenting  a  fairly 
general  computational  scheme  for  a  realistic  reaction  mechanism. 
In  addition  to  the  study  of  the  combustion  process,  a  m.athem.atical 
model  for  the  entire  Otto  cycle  has  also  been  presented.   Tlie  model 
can  be  used  to  study  the  effect  of  different  parameters  on  the  exhaust 
emissions. 


CIIAPTER   II 
SPARK- IGNITION   ENGINE   COMBUSTION  AND   RESULTING   H-IISSIONS 

2.1  Introduction 

In  this  chapter,  the  formation  of  the  various  pollutants,  the 
factors  affecting  them,  and  the  various  control  methods  available  are 
briefly  discussed.   Finally,  a  description  of  the  combustion  process 
together  with  some  introductory  remarks  on  laminar  flames  is  provided. 

TTiere  are  three  important  sources  of  emissions  from  an  auto- 
mobile:  the  engine  exhaust,  the  crankcase  blowby,  and  evaporation 
losses  from  the  fuel  system  vents,  the  engine  exhaust  emitting  by  far 
the  largest  amount  of  pollutants.   Thus  most  of  the  research  on 
automobile  pollution  has  been  concerned  vzith  reducing  the  exhaust 
emissions.  Tne    role  of  the  other  two  sources  of  emissions  is  discussed 
briefly  for  the  sake  of  completeness. 

2.2  Crankcase   Rlov^by 

This  refers  to  the  gases  escaping  from  the  engine  cylinders 
into  the  crankcase  through  the  gaps  between  the  sealing  surfaces  of  the 
piston  and  the  cylinder  v/alls.   This  leakage  which  occurs  mainly  during 
the  expansion  and  the  compression  strokes  increases  with  increasing 
load  (increased  pressures  and  mass  flow)  and  typically  consists  of 
about  85%  unburnt  fuel-air  mixture  and  15%  combustion  products.   In 
terms  of  hydrocarbon  concentrations,  this  is  about  6000-12000  parts 


per  million.   On  uncontrolled  engines  the  crankcase  emissions  account 
for  about  20%  of  the  total  emissions  from  an  automobile. 

Crankcase  emission  control  systems  are  at  present  in  a  well- 
developed  state,  so  that  this  source  of  emissions  has  been  virtually 
eliminated.   All  the  control  systems  use  the  same  basic  approach: 
recycling  of  the  gases  from  the  crankcase  back  into  the  combustion 
chamber,  and  this  is  commonly  referred  to  as  'positive  crankcase 
ventilation'  (PCV) . 

2.3   Evaporation  Losses 

These  losses  occur  from  the  carburetor  and  the  fuel  tank,  and 
on  uncontrolled  engines  account  for  about  18-20%  of  the  total  hydro- 
carbon emissions.   Losses  from  the  carburetor  occur  mainly  after  the 
engine  is  shut  off  when  the  temperature  of  the  fuel  system  rises  due 
to  the  heat  stored  in  the  engine.   The  losses  depend  upon  the  operating 
conditions  imjnediately  before  engine  shutdown,  the  fuel  volatility, 
carburetor  bowl  capacity,  and  the  ambient  temperature. 

llie  losses  from  the  carburetor  are  generally  small  v/hen  the 
engine  is  running,  and  the  fuel  tank  accounts  for  the  majority  of  the 
evaporation  losses.   The  losses  from,  the  fuel  tank  are  dependent  on 
the  fuel  volatility,  its  temperature,  what  fraction  of  the  tank  is 
filled  with  fuel,  and  the  ambient  temperature.   Fuel  evaporation  loss 
control  devices  store  the  vapors  from  the  fuel  tank  and  the  carburetor. 
The  vapors  are  subsequently  burnt  in  the  engine  combustion  chamber.   The 
various  controls  in  use  are  based  on  different  modifications  of  this 
basic  principle.   Evaporation  loss  controls  have  been  in  effect  since 
1972. 


As  mentioned  before,  these  two  losses  (crankcase  and  evaporation) 
contribute  about  20%  each  to  the  total  emissions  of  an  uncontrolled 
vehicle.   The  controls  developed  for  minimizing  the  pollutants  emitted 
thereby  are  in  a  fairly  good  state  of  development  so  that  these  two 
sources  of  emission  have  been  virtually  eliminated.   All  modern  auto- 
mobiles (19  72  and  later  year  models)  are  equipped  with  these  controls 
so  that  presently  the  engine  exhaust  accounts  for  almost  all  the 
pollutants . 

2.4   Formation  of  Pollutants 

We  shall  now  look  at  the  formation  of  the  various  pollutants 
vrhich  occur  in  the  exhaust — our  primary  concern  here. 

2.4.1   Hydrocarbons 

The  presence  of  hydrocarbons  in  the  engine  exhaust  shows  that 
combustion  failed  to  occur  in  certain  parts  of  the  engine  cylinder. 
Calculations  (chemical  equilibrium  or  finite  rate)  reveal  that  the 
concentrations  of  the  hydrocarbons  in  high  temperature  combustion 
products  are  negligibly  sm.all.   The  source  of  the  unburnt  hydro- 
carbons is  not  this  fairly  homogeneous  part  of  the  combustion  chamber. 
The  major  source  of  the  unburnt  hydrocarbons  is  the  quench  layer 
adjacent  to  the  relatively  cool  combustion  chamber  x-7alls ,  wherein  the 
oxidation  of  the  fuel  is  only  partially  completed.   Further,  due  to 
improper  mixing,  the  local  mixture  composition  may  not  be  within  the 
flammability  limits  and  this  would  partly  account  for  the  unburnt 
hydrocarbons.   In  an  Otto  cycle  engine  the  quenching  is  by  far  the 
more  im.portant  of  the  two. 


The  quench  layer  may  be  considered  to  consist  of  two  zones: 
one  knov^rn  as  the  total  quench  zone  which  is  immediately  adjacent  to 
the  wall  and  where  virtually  no  chemical  reactions  occur,  and  the 
other  known  as  the  partial  quench  zone  where  the  temperature  is  high 
enough  for  cool  flame  and  slow  oxidation  reactions  to  take  place,  but 
not  at  a  fast  enough  rate  to  cause  complete  combustion. 

After  flame  propagation,  some  of  the  hydrocarbons  in  the  quench 
layer  mix  with  the  burnt  gases  and  experience  partial  oxidation.   The 
availability  of  oxygen  and  the  degree  of  turbulence  strongly  influence 
the  extent  of  this  oxidation.   Also  all  the  hydrocarbons  present  in 
the  quench  layer  are  not  expelled  during  exhaust  and  some  remain  as 
residuals  in  the  cylinder  for  the  next  cycle.   Reactions  continue  in 
the  gases  e>:pelled  into  the  exhaust  pipe  causing  further  oxidation  of 
the  unburnt  hydrocarbons.   For  these  reasons,  the  amount  of  hydrocarbons 
actually  present  in  the  exhaust  gases  is  considerably  less  than  the 
hydrocarbons  present  in  the  initial  quench  layer. 

Besides  the  quench  layer  the  crevices  between  the  piston  rings 
and  other  dead  spaces  and  pockets  in  the  combustion  chamber  which  prevent 
flam.e  propagation  also  contribute  to  the  exhaust  hydrocarbons. 

The  hydrocarbons  in  the  exhaust  are  a  complex  mixture  of 
oxygenated  and  partially  oxidized  compounds  together  with  products 
produced  by  the  cracking  of  the  fuel;  thus  the  exliaust  hydrocarbons 
are  different  in  composition  from  those  produced  by  crankcase  blox^by 
and  evaporation  losses  which  contain  mainly  the  fuel  itself. 


l.k.l      Carbon  monoxide 

This  is  a  product  of  incomplete  combustion  and  results  when 
sufficient  oxygen  is  not  available  for  combustion,  either  on  an  overall 
or  a  local  basis.   Significant  amounts  of  carbon  monoxide  may  be  present 
in  the  combustion  chamber  immediately  after  the  combustion  is  complete 
even  when  sufficient  oxygen  is  available.   This  can  be  sho^<m  by  a 
chemical  equilibrium  analysis.   However,  as  the  gases  cool  during  the 
expansion  stroke,  much  of  this  carbon  monoxide  gets  converted  to  carbon 
dioxide,  so  that  except  for  fuel-rich  mixtures,  the  amount  of  carbon 
monoxide  in  the  exhaust  may  be  of  the  order  of  about  1%  by  volume  or 
about  10,000  parts  per  million.   The  conversion  of  carbon  monoxide  to 
carbon  dioxide  during  the  expansion  stroke  cannot  be  predicted  from 
chemical  equilibrium  considerations  since  the  chemical  reactions 
involved  are  not  rapid  enough  to  reach  equilibrium.   As  such,  the  actual 
concentration  of  carbon  monoxide  at  exhaust  is  higher  than  that  if 
chemical  equilibrium  v/ere  to  prevail. 

2.4.3   Oxides  of  nitrogen 

A  brief  discussion  is  provided  here;  details  may  be  found  in 
Chapter  V.   The  m.ost  important  of  the  oxides  of  nitrogen  of  concern 
here  is  nitric  oxide  (NO).   The  formation  of  nitric  oxide,  a  high- 
enthalpy  species,  occurs  due  to  the  existence  of  high  temperatures 
during  the  combustion  process.   The  process  of  formation  is  slow 
relative  to  the  time  scale  of  other  processes  occurring  in  a  typical 
internal  combustion  engine,  and  so  has  to  be  considered  on  a  rate 
dependent  basis  and  not  on  a  chemical  equilibrium  basis.   During  the 


expansion  process,  the  concentration  of  nitric  oxide  is  essentially 
frozen  at  an  early  stage  and  hence  the  concentration  at  exhaust  is 
higher  than  that  if  chemical  equilibrium  were  to  prevail. 

2.5   Effects  of  Various  Parameters  on  Engine  Emissions 

2.5.1   Hydrocarbons 

The  hydrocarbon  emissions  depend  primarily  on  the  air-fuel 
ratio  and  engine  load  conditions.   llae  emissions  decrease  with  an 
increase  in  the  air-fuel  ratio  to  a  certain  extent  until  the  mixture 
is  so  lean  that  regular  flame  propagation  is  hindered.   This  is  commonly 
referred  to  as  the  'misfire'  condition.   With  further  increase  of  the 
air-fuel  ratio,  the  hydrocarbon  emissions  increase.   By  increasing  the 
turbulence  and  having  heated  inlet  manifolds,  regular  flame  propagation 
can  be  extended  to  very  lean  mixtures.   In  terms  of  the  operating  condi- 
tions the  influence  of  the  air-fuel  ratio  may  be  summed  up  as:   hydrocarbon 
emissions  are  a  maximum  during  engine  warming,  idle  and  deceleration. 
However,  since  the  amount  of  air  intake  under  these  conditions  is  small, 
the  absolute  amount  (gm./hr.)  of  emJ.ssions  is  not  necessarily  a  maximum. 
The  emissions  also  increase  with  increasing  load  and  com^pression  ratio 
and  early  spark  timing.   With  greater  spark  retard,  combustion  is 
completed  later  during  the  expansion  stroke  and  is  spread  over  a  longer 
time.   The  exhaust  temperatures  are  also  higher,  thus  promoting  oxidation 
in  the  exhaust  system.   The  result  is  lower  emissions,  but  the  performance 
may  be  adversely  affected. 

The  emissions  decrease  with  speed  to  a  certain  extent.   At  very 
high  speeds,  the  ignition  voltage  available  at  the  spark  plugs  decreases 


due  to  the  short  time  available  for  coil  buildup.   This  causes  only 
partial  combustion,  increasing  the  concentrations  of  unburnt  hydro- 
carbons.  The  recently  introduced  capacitive-discharge  ignition 
systems  maintain  a  high  enough  voltage  irrespective  of  the  operating 
speed.   Increasing  the  average  operating  temperature  and  freeing  the 
combustion  chamber  from  deposits  serve  to  decrease  the  hydrocarbon 
emissions. 

The  quenching  effect  and  hence  the  emission  of  hydrocarbons 
are  decreased  by  a  decrease  in  the  surface  area  to  volume  ratio 
within  the  combustion  chamber,  so  that  a  larger  part  of  the  mixture 
is  in  the  flammable  zone.   This  is  achieved  by  increasing  the 
compression  ratio,  increasing  the  stroke  to  bore  ratio  and  m.odifying 
the  geometry  to  obtain  compact  combustion  chambers  avoiding  dead 
spaces  and  pockets.   Increasing  the  combustion  chamber  surface 
temperature  also  serves  to  decrease  the  extent  of  the  quenching  effects, 

2.5.2   Carbon  monoxide 


Like  liydrocarboas ,  the  carbon  monoxide  emissions  also  depend 
on  the  air- fuel  ratio  and  engine  load  conditions.   The  carbon  monoxide 
emissions  are  a  maximum  during  idling  and  deceleration  when  the  air- 
fuel  ratio  is  towards  the  fuel-rich  side.   As  with  hydrocarbons,  the 
absolute  amount  (gm./hr.)  of  the  emissions  under  these  conditions 
need  not  necessarily  be  at  its  maximum  value.   The  carbon  ni.onoxide 
emissions  also  increase  with  load  and  during  acceleration  periods. 
The  concentration  is  generally  at  a  minimum,  during  moderate  power 
cruising.   Other  factors  which  affect  the  hydrocarbon  em.issions 
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also  affect  the  carbon  monoxide  emissions  in  a  similar  manner. 

2.5.3   Oxides  of  nitrogen 

A  detailed  discussion  is  provided  in  Chapter  V;  here  only 
brief  mention  is  made  of  the  effect  of  the  important  parameters. 

The  nitric  oxide  emissions  depend  primarily  on  three  factors: 

(i)     The  emissions  are  a  maximum  at  an  air-fuel  ratio  corresponding 
to  about  10%  excess  air  (as  compared  to  stoichiometric),  and  decrease 
for  richer  or  leaner  mixtures. 

(ii)     Higher  combustion  temperatures  favor  larger  concentrations 
of  nitric  oxide. 

(iii)    llie  emissions  are  greater,  the  greater  the  time  available 

for  the  formation  of  nitric  oxide,  this  being  a  strongly  time -dependent 
process . 

2.6   Emission  Standards 


llie  standards  set  of  19  76  vehicles  are  (based  on  constant 
volume  sample  cold  start  test  average  with  a  constant  volume  sample 
hot  start  test,  both  with  the  Federal  22-minute  driving  cycle): 

0.41  gm. /vehicle-mile  for  hydrocarbons, 

3.40  gm. /vehicle-male  for  carbon  monoxide,   and 

3.00  gm. /vehicle-mile  for  oxides  of  nitrogen. 

Oxides  of  nitrogen  have  to  be  cut  to  0.4  gm. /vehicle-mile  for  19  77 
vehicles . 

For  the  sake  of  comparison,  the  emission  levels  without  any 
exhiaust  controls  for  a  typical  automobile  are  given  below: 
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Hydrocarbons:        850  ppm  (approximately  11  gm. /vehicle-mile) 
Carbon  Monoxide;      3.4%  by  volume  (approximately  80  gm. /vehicle-mile) 
Oxides  of  Nitrogen:   1000  ppm  (approximately  4  gm. /vehicle-mile) 

With  control  techniques  in  current  use  the  emission  levels  are 
approximately  (1973  standards  based  on  constant  volume  sample  cold 
start  test): 

3.4  gm. /vehicle-mile  for  hydrocarbons, 
39.0  gm. /vehicle-mile  for  carbon  monoxide,  and 

3.0  gm. /vehicle-mile  for  oxides  of  nitrogen. 

[Tnese  figures  would  be  somewhat  different  (lower  for  hydrocarbons  and  car- 
bon monoxide)  with  the  testing  procedure  employed  for  the  1976  standards.] 

2.7   Control  Methods 


A  brief  mention  of  the  various  methods  used  to  control  the  exhaust 
emissions  is  made  here. 

(1)  Leaner  air-fuel  ratios  tend  to  decrease  both  the  hydrocarbons  and 
carbon  m.onoxide  emissions;  oxides  of  nitrogen  are  also  reduced  after  a 
certain  air-fuel  ratio.   However,  at  these  air-fuel  ratios  if  regular 
flame  propagation  cannot  occur,  the  former  two  emissions  may  actually 
increase  rather  than  decrease.   Presently,  research  is  being  conducted 
to  study  the  extension  of  the  lean  limit  and  also  how  to  cope  with  the 
generally  adverse  effects  on  driveability  and  some  other  operational 
characteristics  of  the  engine  at  extremely  lean  air-fuel  ratios. 

(2)  Secondary  air  pumped  either  into  the  engine  cylinders  directly 
during  the  expansion  stroke  or  later  in  the  e>diaust  manifold  brings 
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about  a  reduction  in  the  hydrocarbon  and  carbon  monoxide  emissions  by 
promoting  more  nearly  complete  oxidation.  The  nitric  oxide  emissions 
are  not  significantly  affected,  they  may  increase  slightly. 

(3)  Design  modifications:   These  include  improved  carburetor  performance, 
a  fast  acting  and  more  accurate  choke,  an  inductive  or  electronic  igni- 
tion system  with  modified  timing,  reduced  compression  ratios,  air 
preheating,  etc.   Some  of  these  call  for  a  compromise  between  the 
operational  characteristics  and  the  emission  levels.   Previously,  the 
best  design  was  considered  to  be  the  one  that  resulted  in  the  best 
operational  characteristics,  and  this  is  generally  different  from  the 
condition  for  lowest  emission  levels. 

(4)  Post-cylinder  exhaust  treatment  of  the  emissions  prior  to  their 
exit  from  the  exhaust  system  can  be  done  in  several  V7ays  besides 
secondary  air  injection  mentioned  above.   The  most  important  of  these 
is  the  use  of  catalytic  converters  to  promote  the  reactions  in  the 
exliaust  system  in  a  direction  which  favors  the  oxidation  of  the  unburnt 
or  partially  burnt  products,  and  also  to  decompose  nitric  oxide  into 
molecular  nitrogen  and  oxygen,  the  two  processes  generally  being  carried 
out  in  different  stages.   Extensive  research  has  been  done  in  this  area 
but  still  some  problems  remain  before  the  system  can  be  put  in  operation 
on  an  extensive  commercial  basis.   The  high  initial  cost  involved,  the 
Interference  of  the  lead  in  the  gasoline  with  the  catalyst  action,  and 
the  associated  reduction  of  effective  lifetime  are  some  of  the  problems 
to  be  considered. 

Other  modes  of  exhaust  treatment  include  flam.e  afterburners,  but 
these  do  not  result  in  the  reduction  of  nitric  oxide  emissions,  and 
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besides  require  considerable  auxiliary  equipment. 

(5)   Exhaust  gas  recirculation  serves  to  reduce  the  nitric  oxide  emissions, 
but  the  carbon  monoxide  and  hydrocarbon  emissions  are  adversely  affected 
as  are  some  of  the  operational  characteristics  of  the  engine. 

There  are  several  other  control  techniques  which  have  been 
suggested  and  the  reader  is  referred  to  the  appropriate  literature 
on  the  subject. 

2.8  Description  of  the  Combustion  Process  in  an  Otto  Cycle  Engine 

In  a  spark  ignition  engine,  the  high  temperature  of  the  electric 
spark  causes  the  ignition  of  a  small  portion  of  the  fuel-air  mixture 
in  the  immiediate  vicinity  of  the  spark  plug  gap.   The  result  is  a  m.ore 
or  less  spherical  volume  of  hot  gas  consisting  of  burnt  products  which 
are  rapidly  equilibrated  due  to  the  reactions  being  very  fast  at  the 
high  temperature  involved.   Tlais  causes  a  nearly  spherical  flame  to 
originate  from  the  point  of  ignition.   Spherical  flames  are  truly  one- 
dim.ensional  systems  by  virtue  of  their  geometry,  whereas  flat  flames  are 
quasi-one-dimensional  systems,  there  being  small  but  negligible  changes  in 
the  various  state  properties  in  the  direction  perpendicular  to  that  of 
propagation.   Ihe  fuel-air  mixture  in  the  com.bustion  chamber  is  exposed 
to  the  flame,  layer  by  layer,  until  combustion  is  com.plete  as  far  as 
practicable.   In  practice,  combustion  is  never  totally  complete  due  to  the 
quenching  effects  of  the  walls  and  other  reasons.   Tne  entire  combustion 
process  is  extremely  rapid,  e.g.,  for  combustion  spread  over  40  degrees 
of  crank  revolution  at  2000  revolutions  per  minute,  the  process  x^7ould  be 
completed  in  approxim.ately  3.33  milliseconds.   In  spite  of  the  rapidity  of  the 
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process,  at  any  given  Instant  of  time,  the  fraction  of  the  total 
fuel-air  mixture  actually  being  burnt  is  not  significant  if  normal 
combustion  is  occurring.   Tliis  would  not  be  the  case  when  detonation 
and  knocking  occur.   Therefore,  under  normal  circumstances,  the 
combustion  process  occurs  with  uniform  release  of  energy  and  may  be 
approximated  in  the  initial  stages  by  the  model  of  a  one-dimensional 
laminar  spherical  flame  propagating  into  a  premixed  fuel-air  mixture. 
The  assumption  of  a  laminar  flame  is  made  to  render  the  problem 
solvable  (in  general,  turbulence  would  be  present  in  varying  degrees), 
there  being  no  convenient  numerical  techniques  to  handle  a  turbulent 
flame  for  a  general  case.   As  the  flame  propagates,  it  would  deviate 
from  its  spherical  shape,  the  deviation  depending  on  the  geometry  of 
the  combustion  chamber  walls  with  respect  to  the  point  of  ignition. 
This  deviation  can  be  accounted  for  by  considering  the  general  three- 
dimensional  flame  equations  in  their  time-dependent  form,  and  a  solution 
so  obtained  would  give  the  flame  shape  and  position  as  a  function  of 
time.   Even  with  modem  electronic  computers  this  problem  is  almost 
insurmountable,  and  would  hardly  give  more  insight  to  the  actual 
combustion  process  than  the  solution  of  a  one -dimensional  problem. 

2.9   Laminar  Flames:   Introductory  Remarks 

Laminar  premixed  flames  represent  a  unique  case  of  vzave 
propagation — the  propagation  of  a  self-sustaining  exothermic  reaction 
wave  into  a  mixture  of  combustible  gases.   Heat  conduction,  diffusion, 
and  a  set  of  rapid  exothermic  chemical  reactions  all  are  significant 
in  the  propagation  of  such  a  flame.   Thus  this  is  one  of  the  earliest 
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combustion  problems  that  has  necessitated  the  simultaneous  consideration 
of  both  the  fluid  dynamic  and  chemical  kinetics  aspects  for  its  solution. 
Laminar  flame  velocity  (also  referred  to  variously  as  the  burning 
velocity  or  the  normal  combustion  velocity),  which  is  of  primary 
importance  in  the  study  of  premixed  flames,  is  defined  as  the  velocity 
at  which  unburnt  gases  move  through  a  combustion  wave  in  the  direction 
normal  to  the  wave  surface.   This  definition  refers  to  the  velocity 
of  unburnt  gas  relative  to  the  flame  front  as  determined  at  some  point 
far  enough  removed  from  the  flame,  such  that  the  influence  of  the  flame 
itself  at  the  point  is  negligible. 

The  propagation  characteristics  of  steady,  adiabatic,  one- 
dimensional  laminar  premixed  flames  are  governed  primarily  by  two 
properties  of  the  system — the  activation  energy  of  the  rate  controlling 
step  and  the  Lewis  number,  where  the  Lewis  number  represents  the  ratio 
of  the  energy  transported  by  conduction  to  that  transported  by  diffusion. 
Tne  geometry  is  of  secondary  importance,  and  therefore,  within  reason- 
able limits,  the  laminar  flame  velocity  may  be  defined  independent  of 
a  particular  experimental  apparatus.   It  will  be  seen  later  that  the 
flame  velocity  is  an  eigenvalue  of  the  conservai.j-on  equations. 

A  laminar  flame  is  made  up  of  a  series  of  fairly  well-defined 
regions.   The  thickness  of  the  flame  zone  is  rather  arbitrary,  similar 
to  the  thickness  of  a  boundary  layer.   In  the  first  region,  the  unburnt 
gas  undergoes  appreciable  changes  in  composition  and  temperature  due  to 
the  transport  processes  of  diffusion  and  thermal  conduction,  but  the 
reaction  rates  are  negligible.   This  is  referred  to  as  the  preheat 
zone.   After  this  there  is  a  relatively  thin  primary  reaction  zone 
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wherein  all  the  important  rapid  chemical  reactions  occur.   This  is  the 
luminous  zone  in  hydrocarbon  flames.   In  addition  to  this,  there  may 
be  one  or  more  secondary  reaction  zones  vjhose  spatial  dimensions  are 
an  order  of  magnitude  larger  than  those  of  the  primary  reaction  zone, 
and  wherein  slow  reactions  (e.g.,  carbon  monoxide  afterburning  in 
hydrocarbon  flames)  and  radical  recombination  reactions  occur.   At 
low  enough  pressures  the  radical  recombination  reactions  and  the 
afterburning  may  occur  in  distinctly  separate  zones;  on  the  other  hand, 
with  increasing  pressure,  the  secondary  reaction  zone  may  partly  merge 
into  the  primary  zone.   The  existence  of  these  zones  and  carbon  monoxide 
afterburning  has  been  confirmed  by  experiments. 

2.10   Closing  Pvemarks 

An  introduction  to  the  study  of  the  combustion  process  in  and 
the  emissions  from  an  Otto  cycle  engine  has  been  presented.   These 
ideas  will  be  developed  further  in  the  succeeding  chapters. 

Chapter  III  presents  a  review  of  the  extensive  literature  on 
the  methods  of  solution  of  the  laminar  flame  propagation  equations. 

Chapter  IV  presents  a  mathematical  analysis  of  the  combustion 
process — the  form.ulation  and  solution  of  the  laminar  flame  propagation 
equations . 

In  Chapter  V  a  mathematical  model  has  been  developed  to  obtain 
a  detailed  thermodynamic  analysis  of  the  Otto  cycle. 

Results  and  discussions  are  presented  in  Chapter  VI  and  finally 
Chapter  VII  is  used  for  some  closing  remarks. 


CHAPTER  III 
LITERATURE  REVIEW 

The  problem  of  laminar  flame  propagation  in  premixed  gases  has 
been  studied  fairly  extensively.   The  problem  was  first  studied  by 
Mallard  and  LeChatelier  (1883).   Their  theory  has  been  referred  to  as 
the  thermal  theory  in  the  sense  that  they  considered  thermal  effects  to 
be  of  primary  importance  and  the  chemical  kinetics  to  be  secondary. 
Later  on,  Taffanel  (1913)  and  Daniell  (19  30)  demonstrated  for  a  highly 
simplified  model  that  the  flame  propagation  velocity  is  proportional 
to  the  square  root  of  the  ratio  of  the  thermal  conductivity  to  the 
specific  heat  at  constant  pressure  and  to  the  square  root  of  the  rate 
of  the  reaction. 

An  excellent  review  of  the  early  literature  (up  to  about  1950) 
on  the  steady  state  flame  propagation  in  premixed  gases  has  been 
presented  by  Evans  [1].   Williams  [2,  p.  95]  has  presented  a  review 
with  references  to  more  recent  literature  (up  to  about  1950). 

The  assumption  of  Lewis  number  equal  to  unity  (normal  diffusion) 
has  been  widely  made  in  flame  analysis.   This  assum.ption  implies  that 
the  energy  fluxes  due  to  diffusion  and  conduction  balance  each  other, 
and  so  the  enthalpy  is  constant  throughout  the  flame  zone;  hence  the 
enthalpy  equation  need  not  be  integrated.   lliis  simplifies  the  problem 
considerably  and  even  permits  closed  form  solutions  for  a  few  simple 
cases.   The  assumption  of  unity  Lewis  number  is,  however,  not  valid  for 
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many  gases  encountered  in  combustion  problems,  two  notable  examples 
being  hydrogen-air  and  propane-air  mixtures  (Le  =  3.30,  and  Le  =  0.49, 
respectively) . 

Assuming  that  the  activation  energy  of  the  rate  controlling 
step  and  the  heat  of  reaction  are  large  and  that  all  the  reaction  occurs 
at  the  hot  boundary,  Zeldovich  and  Frank-Kamenetski  [3]  were  able  to 
obtain  a  closed  form  solution  for  the  burning  rate  eigenvalue.   Semenov 
[Sokolik,  4]  has  modified  their  solution  for  non-normal  diffusion  by 
multiplying  the  flame  velocity  as  calculated  by  the  Zeldovich-Frank- 
Kamenetski  formulation  by  the  factor  1/v^^e  for  first  order  reactions 
and  by  1/Le  for  second  order  reactions.   This  analysis  is  also  limited 
by  the  requirement  of  large  or  infinite  activation  energies.   Von  Karman's 
zeroth  and  first  order  approximations  [5]  also  yield  closed  form  solu- 
tions, the  former  being  a  special  case  of  the  Zeldovich-Frank-Kamenetski 
analysis  wherein  only  the  first  term  of  a  series  expansion  for  the 
eigenvalue  is  retained.   Some  of  the  other  analytical  solutions  are 
those  of  Boys  and  Corner  [6,  7]  and  Adams  [8].   Wilde  [9]  empirically 
modified  Adams'  solution  thereby  obtaining  better  agreement  with  'exact' 
numerical  solutions,   Hirschfelder  [10]  has  obtained  two  approximate 
solutions  in  one  of  which  the  eigenvalue  is  obtained  as  a  root  of  a 
quadratic  equation  and  in  the  other  as  a  root  of  a  third  order  equation. 
Tlie  accuracy  is  comparable  to  that  obtained  with  Wilde's  solution. 
Williams  [2]  has  given  a  graphical  comparison  of  the  accuracy  of  these 
approximate  methods. 

Spalding  [11]  has  proposed  a  'centroid  rule'  v/hich  provides  a 
convenient  way  of  estimating  the  flame  speed  to  an  accuracy  of  about  1% 
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for  the.  case  in  which  the  reaction  rate  is  an  explicit  function  of 
temperature.   This  is  true  for  unity  Lewis  number  of  the  reactants 
and  radicals  for  the  following  class  of  reactions:   (a)  single  step 
reactions,  (b)  reactions  satisfying  the  steady  state  approximation, 
and  (c)  explosive  chain  reactions.   The  method  gives  an  expression 
for  the  flame  speed  in  terms  of  the  centroid  of  the  reaction  rate 
function  without  resorting  to  the  solution  of  the  governing  dif- 
ferential equations.   Graphs  are  provided  for  estimating  the  flame 
speed  for  non-unity  Lewis  number.   These  show  that  when  the  Lewis 
number  of  the  radicals  is  greater  than  unity  the  flame  speed  in- 
creases, but  only  slightly.   Small  values  of  Lewis  number  (of  the 
radicals),  on  the  other  hand,  greatly  reduce  the  flame  speed.   For 
the  reactants,  a  Lewis  number  greater  than  unity  decreases  the 
flame  speed. 

Rosen  [12,  13]  has  suggested  a  variational  method  employing 
the  Rayleigh-Ritz  method  for  obtaining  successively  better  estimates 
of  the  functional  and  eigenvalue.   He  has  also  given  a  plausible 
derivation  of  the  centroid  rule  proposed  by  Spalding  [11].   The 
centroid  rule  method  has  been  extended  for  the  case  of  non-normal 
diffusion  with  bimolecular  reactions  (A+B-^C)  in  a  stoichiometric 
two-component  system  by  Adler  [14].   The  flame  speed  has  been  sho\ra  to 
be  inversely  proportional  to  the  Lewis  number  of  the  reactants  for 
large  activation  energies. 

Spalding  [15]  has  compared  several  approximate  solutions  for 
the  case  of  unity  Lewis  number  and  temperature  explicit  reactions  and 
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arranged  them  in  order  of  merit  as:   Wilde  (best),  Adams,  Von  Karman- 
Penner,  Zeldovich-Frank-Kamenetski,  and  Corner.   Wilde's  method  has 
been  extended  to  study  the  effect  of  the  diffusion  coefficient  of 
reactants  and  radicals  on  the  flame  speed  and  fairly  good  agreement 
with  'exact'  solutions  is  reported. 

Adler  [16]  has  considered  temperature  explicit  reaction  rate 
functions  and  obtained  limits  on  the  burning  rate  eigenvalue  independent 
of  the  type  of  reaction  considered  for  a  fixed  centroid  of  the  reaction 
rate  function.   For  the  limiting  case  of  high  activation  energy  the 
upper  and  lov7er  limits  are  close,  so  that  the  flame  speed  may  be 
obtained  with  good  accuracy. 

Favin  and  Westenberg  [17]  have  considered  one-dimensional 
steady  premixed  spherical  laminar  flames  and  solved  the  equations  for 
the  case  of  unity  Lewis  num.ber  and  a  unimolecular  reaction  process. 
They  have  integrated  the  equations  from  the  cold  to  the  hot  boundary 
and  considered  the  flam.eholder  temperature,  rather  than  the  m.ass  flow 
rate,  as  the  eigenvalue. 

Adler  and  Spalding  [18]  have  obtained  analytical  solutions  for 
the  temperature  gradient,  radical  concentrations,  and  the  flame  speed 
for  the  case  of  non-branching  chain  reactions  with  negligible  chain 
breaking. 

Ue  Sendagorta  [19]  has  obtained  a  numerically  exact  solution 
for  one  value  of  the  reduced  activation  temperature  E/RT^  and  used 
this  to  shoxs?  that  good  accuracy  could  be  obtained  by  using  an 
exponential  approxim.ation  for  the  species  flux  fraction  in  analytical 
solutions.   His  solution  covers  both  first  and  second  order  single 
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step  reactions  and  is  considerably  more  accurate  than  the  other 
analytical  solutions  cited  above.   Non-normal  diffusion  has  been 
considered  by  Wilde  [9],  de  Sendagorta  [19]  and  Von  Karraan-Penner 
[20]  utilizing  Karman-Pohlhausen  profile  methods. 

Millan  and  da  Riva  [21]  have  obtained  a  solution  for  the 
distribution  of  radicals  in  laminar  flames  based  on  the  steady  state 
approximation.   A  general  criterion  for  the  applicability  of  this 
approximation  is  given  and  it  is  shown  that  this  assumption  is  a  first 
order  approximation  to  the  actual  solution,  and  this  solution  can  be 
improved  by  a  method  of  successive  approximations. 

Spalding  and  Adler  [22,  23]  have  considered  one-dim.ensional 
laminar  flames  with  an  enthalpy  gradient  (defined  as  positive  for 
heat  loss  at  the  cold  boundary)  which  is  useful  when  heat  losses  by 
conduction  or  radiation  are  significant,  and  obtained  numerical 
solutions.   For  unity  Lewis  number,  they  have  shown  that  the  flame 
speed  decreases  or  increases  as  the  enthalpy  gradient  decreases  or 
increases  from  zero,  the  effect  being  more,  the  higher  the  activation 
energy.   The  flame  speed  decreases  to  zero  for  a  certain  critical 
negative  value  of  the  enthalpy  gradient. 

Heat  losses  and  flame  propagation  limits  due  to  heat  losses 
have  also  been  considered  in  the  analyses  presented  by  Mayer  [24], 
Spalding  [25],  Adler  [26],  and  /dler  and  Kennerley  [27].   From  these 
analyses  it  can  be  concluded  that  for  a  given  amount  of  heat  loss, 
a  given  combustible  mixture  has  two  possible  flam.e  speeds,  the  lower 
one  of  which  is  normally  unstable.   With  increasing  heat  loss,  the 
two  speeds  approach  each  other  and  finally  coincide  for  a  critical 
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value  of  the  heat  loss.   For  still  greater  heat  loss,  the  flame  speed 
is  imaginary.   The  critical  value  of  heat  loss  is  identified  with  the 
flammability  limit.   Some  other  studies  of  non-adiabatic  flames  are 
given  in  References  28,  29,  and  30. 

With  the  advent  of  modern  electronic  computing  machines,  the 
emphasis  has  gradually  shifted  to  obtain  more  meaningful  solutions  with 
less  drastic  assumptions,  by  using  numerical  methods.   These  studies 
involve  less  stringent  assumptions,  are  more  accurate,  and  include 
better  estimates  of  thermodynamic  and  transport  properties,  and  also 
consider  the  effects  of  chain  reactions  instead  of  the  classical 
unimolecular  process  R  ^  P  (R  =  reactant,  P  =  product),  considered 
in  most  of  the  previous  analytical  solutions.   Firstly,  brief  mention 
is  made  of  the  early  numerical  solutions. 

Von  Karman's  second  approximation  [5]  involves  an  iterative 
procedure  with  numerical  integrations.   The  approximate  solution  of 
Boys  and  Corner  [6]  can  be  improved  by  numerical  iterations.   Some 
of  the  other  iterative  procedures  requiring  numerical  integrations 
are  those  due  to  Klein  [31]  and  Johnson  and  Nachbar  [32].   The  latter 
one  is  the  most  accurate  of  the  simple  approximate  procedures  and  it 
permits  the  determination  of  the  upper  and  lower  bounds  for  the 
eigenvalue.   Using  their  iteration  scheme  successive  upper  and  lower 
bounds  can  be  obtained  until  the  eigenvalue  converges  to  any  desired 
accuracy.   Generally  the  first  approximation  is  accurate  enough  for 
practical  purposes  so  that  it  is  seldom  necessary  to  use  the  iterative 
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Friedman  and  Burke  [33]  have  obtained  numerical  solutions  for 
the  case  of  a  hypothetical  first  order  irreversible  reaction  with 
exponential  dependence  of  rate  on  temperature  over  a  wide  range  of  two 
param.eters:   dimensionless  activation  temperature  (E/RT  )  and  the 
Lewis  number. 

The  theoretical  analysis  presented  by  Giddings  and  Hirschfelder 
[ 3A ]  was  the  first  one  to  consider  branched  chain  reactions.   The 
assumption  of  a  steady  state  condition  was  made  to  approximate  the 
radical  concentrations  and  the  governing  equations  were  solved  by  the 
integral  equation  approach  devised  by  Klein  [31]. 

Menkes  [35]  has  solved  the  problem  for  a  cylindrical  flame 
by  transforming  the  differential  equation  into  an  integral  equation 
which  is  solved  numerically  by  a  method  of  successive  approxim.ations . 
Tlie  cylindrical  flame  was  investigated  to  get  an  insight  into  the 
stability  analysis.   It  was  found  that  the  flame  radius  decreases 
and  the  specific  mass  flow  increases  rapidly  v;ith  increasing  inlet 
temperature  and  at  constant  flame  thickness.   However,  if  the  flame 
radius  is  decreased,  keeping  the  inlet  temperature  constant,  the 
flame  thickness  increases  rapidly  resulting  in  a  decrease  of  the 
specific  mass  flow  rate.   Similar  results  have  been  obtained  by 
Menkes  [36]  by  an  approximate  analytical  solution.   Jain  and  Ebenezer 
[37]  have  obtained  analytical  solutions  for  a  cylindrical  f lam.e ,  wherein 
the  reaction  rate  is  represented  by  a  step  function,  and  have  studied 
the  accuracy  of  the  thin-flame"  approxim.ation. 


A  thin  flam.e  is  one  whose  thickness  is  small  compared  to  its 
radius  of  curvature,  so  that  the  mass  rate  of  burning  per  unit  area 
can  be  taken  as  equal  to  that  in  a  plane  flame. 
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Spalding  [38]  has  obtained  numerical  solutions  for  the  case 
of  a  steady  laminar  spherical  flame  by  stepwise  integration  from  the 
cold  boundary.   The  model  used  consists  of  a  combustible  gas  emerging 
from  a  point  source  into  an  atmosphere  of  adiabatic  combustion  products. 
The  results  have  been  compared  with  those  obtained  from  the  thin-flame 
theory  and  the  two  are  shown  to  be  asymptotically  equal  for  large  mass 
flow  rates  and  high  reaction  rates.   For  the  same  model  Spalding  and 
Jain  [39]  have  obtained  analytical  solutions  for  the  cases  of  (a)  step 
function  reaction  rate  curves  and  (b )  Adams-type  reaction  rate  curves 
based  on  the  assumption  of  a  constant  temperature  gradient  in  the 
reaction  zone.   It  is  shown  that  the  radius  of  the  flame  depends  more 
on  the  position  of  the  centroid  of  the  reaction  rate  curve  rather 
than  on  its  specific  shape.   The  same  results  have  been  confirmed  by 
m.eans  of  a  resistance  network  analog  comiputer  by  Spalding,  Jain  and 
Samain  [ 40 ] . 

Jain  and  Kumar  [41]  have  employed  the  method  of  matched 
asym.ptotic  expansions  and  obtained  analytical  solutions  for  the 
limiting  case  of  large  activation  temperatures.   They  have  also 
obtained  numerical  solutions  for  a  wide  range  of  values  of  both  the 
reduced  activation  temperature  and  Lewis  number,  for  both  first  and 
second  order  reactions.   Based  on  the  numerical  solutions,  a  reason- 
ably accurate  rule  for  the  flame  speed  eigenvalue  in  terms  of  Lewis 
number  and  the  centroid  of  the  reaction  rate  function  is  suggested; 
it  is  linear  in  Lewis  number  for  first  order  reactions  and  quadratic 
in  Lewis  number  for  second  order  reactions.   There  is  good  agreement 
of  these  results  with  those  of  Spalding  [11,  15]  and  de  Sendagorta  [19], 
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However,  the  present  solution  agrees  with  Adler's  [14]  solution  only 
for  Lewis  numbers  around  unity,  and  there  is  considerable  disagreement 
when  the  Lewis  number  differs  significantly  from  unity  for  which  a 
possible  explanation  has  been  given. 

Besides  the  various  methods  mentioned  above,  other  recent 
methods  have  sought  not  just  to  determine  the  laminar  flame  speed 
but  also  the  detailed  flame  structure  regarding  the  distribution  of 
temperature  and  the  various  chemical  species  over  the  thickness  of 
the  flame  zone.   A  detailed  analysis  of  the  flame  structure  gives 
a  better  insight  into  the  mechanism  of  flame  propagation  and  also 
permits  the  study  of  chemical  reactions  from  the  viewpoint  of 
mechanism  and  kinetics.   The  difficulties  involved  in  such  a  detailed 
analysis  of  the  flame  zone  may  be  summed  up  into  three  categories: 
(a)  a  reliable  postulation  of  the  reaction  mechanism  for  the  oxidation 
of  the  fuel;  (b )  a  reliable  estimate  of  the  transport  properties, 
especially  at  temperatures  encountered  in  flames  and  for  free 
radicals;  and  (c)  the  absence  of  suitable  numerical  techniques  to 
handle  the  complicated  system  of  equations  in  a  convenient  way. 

Basically  two  mathematical  formulations  have  been  widely  used 
in  the  detailed  analysis  of  the  flame  zone.   The  first  one  which  is 
a  trial  and  error  method  is  due  to  Hirschfelder  [42].   This  method 
consists  of  setting  up  the  differential  equations  for  the  conservation 
of  mass  and  energy  and  the  diffusion  equations  in  a  coordinate  system 
which  is  at  rest  with  respect  to  the  flame  front.   The  boundary  conditions 
are  specified  at  the  burnt  and  unbumt  ends.   In  practice,  integrating 
these  equations  to  obtain  the  concentrations  of  species  and  state 
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properties  in  the  flame  zone  requires  a  knowledge  of  the  flame  speed. 
Since  this  is  not  known  a  priori,  the  problem  reduces  to  an  eigenvalue 
problem;  a  solution  to  the  equations  exists  for  only  one  unique  value 
of  the  flame  speed,  known  as  the  eigenvalue  of  the  conservation 
equations  or  the  eigenvalue  of  the  two-point  boundary  value  problem. 
Assuming  various  flame  speeds,  the  equations  are  numerically  integrated 
from  the  burnt  to  the  unburnt  end.   The  solution  which  agrees  best  with 
the  boundary  conditions  is  the  one  assumed  to  be  correct.   Convergence 
criteria  for  this  method  are  not  available.   For  more  complicated  cases, 
where  the  chemical  kinetics  is  governed  by  a  set  of  reaction  steps,  the 
assumption  of  flame  speed  alone  is  not  sufficient  to  start  the  integra- 
tion of  the  conservation  equations ,  and  one  is  confronted  with  a  multi- 
eigenvalue  problem  in  non-separable  equations  [1,  43,  44].   The  number 
of  eigenvalues  is  equal  to  the  number  of  linearly  independent  species 
flux  fractions  in  the  flame  equations.   This  multi-eigenvalue  problem 
has,  however,  been  reduced  to  the  determination  of  a  single  eigenvalue, 
without  loss  of  generality,  in  a  method  introduced  by  Campbell,  Heinen 
and  Schalit  [44].   Tliey  have  obtained  numerical  solutions  for  a  one 
dimensional  flame  with  idealized  kinetics  and  transport  properties. 

Tnree  examples  using  this  trial  and  error  method  have  been  given 
by  Hirschfelder  [42]:   (a)   the  unimolecular  decomposition  of  hydrazine, 
(b)  the  bimolecular  decomposition  of  nitric  oxide,  and  (c)  the  two  step 
chain  m.echanism  describing  the  decomposition  of  ozone.   Agreement  with 
extierimsntal  results  is  fair''"  for  the  first  and  third  cases;  no 


"it  should  be  noted  that  the  solution  method  itself,  even  if 
extremely  powerful,  would  not  give  results  in  good  agreement  with  the 
experimental  data  if  the  accuracy  of  the  reaction  kinetics  and  transport 
property  data  is  poor.   Thus  the  judgment  of  accuracy  should  rather  be 
based  on  the  solution  of  the  same  problem  with  the  same  input  data  by 
different  methods  and  comparison  with  experiments. 
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experimental  data  wei-e  available  for  comparison  in  the  second  case. 

The  second  method  is  due  to  Spalding  [45]  and  consists  of 
setting  up  the  differential  equations  for  conservation  of  mass  and 
energy  and  the  diffusion  equations  in  their  time-dependent  or  unsteady 
form.   Assuming  arbitrary  initial  profiles  for  the  state  properties 
and  the  species  concentrations,  the  equations  are  solved  by  forward 
marching  finite  difference  methods  until  the  various  profiles  of  the 
species  concentrations  and  the  state  properties  converge  to  a  steady 
state.   A  big  advantage  of  this  method  is  that  almost  any  arbitrary 
initial  profile  can  be  used  for  starting  the  solution;  however,  by 
using  profiles  which  closely  resemble  their  steady  state  shape,  the 
computation  can  be  cut  down  considerably.   For  the  temperature  and 
the  major  species,  S-shaped  profiles  asymptotically  approaching  the 
boundary  values   are  com.roonly  used.   The  radicals  are  usually  assumed 
to  be  at  zero  concentration  initially.   The  assumed  profiles  will 
develop  into  the  steady  case  if  enough  energy  is  stored  in  them. 
If  this  initial  energy  is  insufficient  then  the  gradients  of  tem.per- 
ature  and  species  concentrations  will  decrease  m.onotonically  and  the 
flame  is  extinguished.   Thus  this  method  can  also  be  used  to  obtain 
an  estimiation  of  the  minimum  ignition  energy,  if  required.   No 
iterations  are  required  with  this  approach.   The  computational 
advantage  or  disadvantage  of  either  method  has  not  been  conclusively 
established  as  yet. 

This  m.ethod  has  been  demonstrated  by  Spalding  for  the  case  of 
the  hydrazine  decomposition  flame.   Tlie  conservation  equations  reduced  to 
two  partial  differential  equations  and  these  were  solved  sim.ultaneously 
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by  Schmidt's  graphical  method  for  the  solution  of  the  unsteady  heat 
conduction  equation.   Calculations  were    carried  out  for  two  different 
flame  temperatures  and  fairly  good  agreement  between  theoretical  and 
experimental  values  were  obtained. 

Spalding's  method  has  been  used  by  various  researchers  for 
solving  different  flames.   The  later  techniques  used  are  basically 
refinements  of  the  finite  difference  solutions  and  the  inclusion  of 
chain-branching  reactions  and  other  improvements.   A  review  of  some 
of  the  recent  work  using  this  method  is  presented  below. 

Zeldovich  and  Barenblatt  [46]  have  reviewed  previous  work 
on  the  steady  state  propagation  of  flat  flames.   They  have  also 
considered  the  effects  of  the  Lewis  number,  heat  loss  and  chain 
reactions  on  the  velocity  of  propagation  of  flam.es  using  Spalding's 
time-dependent  approach.   Using  a  perturbation  technique  the  stability 
of  the  steady  state  solution  has  been  examined.   They  have  shown  that 
for  non-unity  Lev7is  number,  the  total  enthalpy  reaches  a  maximum  or 
a  minimum  in  the  flame  zone,  depending  on  whether  the  Lewis  number  is 
greater  or  smaller  than  unity.   The  same  result  was  also  obtained 
earlier  by  Lewis  and  Von  Elbe  [47].   Heat  loss  was  foimd  to  decrease 
the  velocity  of  propagation  and  there  was  found  to  exist  a  certain 
critical  value  of  the  heat  loss  parameter  'b '  [q  =  b(t-t^)],  above 
which  a  steady  state  solution  did  not  exist. 

Adams  and  Cook  [48]  have  tised  the  tim.e-dependent  m.ethod  to 
study  the  hydrazine  decomposition  flame  which  has  been  fairly 
extensively  studied.   As  before,  the  problem  V7as  reduced  to  the 
simultaneous  solution  of  two  partial  differential  equations:  one 
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involving  the  temperature  and  the  other  the  mass  fraction  of  hydrazine. 
A  central  difference  formula  was  used  for  approximating  the  space 
derivatives  and  a  f onward  difference  formula  for  the  time  derivatives. 
The  equations  were  then  numerically  solved  by  explicit  forward  marching 
steps  in  time  until  the  profiles  of  temperature  and  mass  fraction  reached 
a  steady  state.   Considering  chain  breaking  reactions  of  both  second  and 
third  order,  it  was  shown  that  the  dependence  of  flame  speed  on  temper- 
ature is  insensitive  to  changes  in  pressure  and  the  mechanism  of  the 
chain  breaking  reactions,  while  the  dependence  on  pressure  is  sensitive 
to  the  mechanism  assumed. 

Dixon-Lewis  [49]  has  implemented  Spalding's  time-dependent 
method  for  the  hydrogen-ox^'gen  flame.   To  reduce  the  computational 
work,  some  simplifications  regarding  the  transport  and  therm.odynamic 
properties  were  made.   The  chemical  kinetics  was  treated  in  terms  of 
a  set  of  elementary  reaction  steps  and  thus  this  analysis  avoids  some 
of  the  misleading  results  of  previous  analyses  which  considered  the 
chemi.cal  kinetics  to  be  governed  by  a  single  reaction  or  a  simple 
unrealistic  mechanism.   Important  conclusions  regarding  the  influence 
of  some  parameters  on  the  hydrogen-oxygen  flame  have  been  made  in 

this  analysis.   For  the  sake  of  computational  stability,  the  standard 

1 
2 


2   1 
stability  condition  of  At/(A!]j)  <   -y  for  the  heat  conduction  equation 


2 
had  to  be  modified  to  At/ (A'J;)   ^  1/30  for  this  problem. 

Detailed  transport  property  calculations  have  been  made  [50] 

and  considered  in  the  analysis  of  the  hydrogen-oxygen-nitrogen  flame 

by  Dixon-Lewis  [51]  using  the  Runge-Kutta-lIerson  procedure  (Runge-Kutta 

method  '.s?ith  variable  step  size). 
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Spalding  et  al.  [52]  have  utilized  the  mathematical  similarity 
of  the  equations  governing  steady  two-dimensional  boundary  layer  flows 
with  the  one-dimensional  unsteady  flame  equations.   The  computational 
procedure  developed  by  Patankar  and  Spalding  [53]  for  the  study  of  the 
hydrodynamical  and  heat  transfer  aspects  of  steady  two-dimensional 
boundary  layer  flows  has  been  adapted  to  the  study  of  the  one-dimensional 
unsteady  hydrazine  decomposition  flame.   Both  these  phenomena  are  governed 
by  sets  of  simultaneous  parabolic  partial  differential  equations.   Tlie 
assumptions  of  unity  Lewis  number  and  constant  and  uniform  specific 
heats  were  made  in  the  solution  method,  although  neither  was  necessary 

for  the  computational  procedure.   Initial  concentration  profiles  of  the 

3      4     5 
form  y  =  y    (lOw  -  15w  +  6w  )  were  used  for  the  major  species,  and 

the  radicals  were  assumed  to  be  initially  at  zero  concentration.   Here, 
w  =  (.i'-''<P,)/(ii     -  ip,)  ,   where  ij;  =  stream  function,  and  the  subscripts  b 
and  u  refer  to  the  burnt  and  unburnt  ends,  respectively.   A  variable 
step  size  for  the  space  dimension  was  used  and  the  equations  were  solved 
by  an  implicit  method  after  linearization  at  each  step.   In  the  absence 
of  formal  stability  criteria,  the  largest  time  step  size  was  found  by 
trial  and  error. 

Wilde  [54]  has  treated  this  problem  as  a  boundary  value  problem 
and  formulated  it  in  three  different  ways:   (a)  a  set  of  first  order 
time- independent  ordinary  differential  equations  for  species  conserva- 
tion, diffusion  and  thermal  flux;  (b)  a  sat  of  second  order  time- 
independent  ordinary  differential  equations;  and  (c)  a  set  of  time- 
dependent  parabolic  partial  differential  equations.   The  first  for- 
mulation was  solved  by  quasi-linearization,  also  known  as  Newton- 
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Raphson-Kantorovich  linearization,  and  then  integrating  by  using 
algorithms  suitable  for  stiff  equations.   The  second  formulation 
was  also  linearized  and  then  the  equations  were  written  in  finite 
difference  form,  the  result  being  a  set  of  linear  algebraic  equation, 
block  tridiagonal  in  form.   For  the  last  formulation,  explicit  finite 
difference  methods  x^7ere  found  to  require  a  very  small  time  step  size 
due  to  the  presence  of  the  term  involving  the  rate  of  production  of 
species  by  chemical  reactions.   Considerable  time  saving  was  obtained 
by  using  forward  difference  formulas  in  time,  implicit  representation 
of  the  space  derivatives  and  an  averaging  of  the  chemical  rate  term 
at  the  present  and  future  tim.e  step.   The  second  formulation  was 
faster  per  iteration  than  the  first  one,  but  required  closer  initial 
estimates  and  a  greater  number  of  iterations  for  convergence.   The 
last  method  was  longer  in  running  time  than  the  other  two  but  was 
surer  in  convergence.   Three  flam.es  were  studied  using  this  method: 
ozone  decoraposition,  hydrogen-bromine  flame,  and  the  hydrogen-oxygen 
flame.   The  calculated  flame  velocities  uere   about  tv/ice  the  cor- 
responding experimental  values.   The  discrepancy  was  explained  by 
blaming  the  unreliable  kinetic  and  transport  data. 

Bledjian  [55]  has  provided  a  general  formulation  which  allows 
for  the  consideration  of  spherical  and  cylindrical  flames  in  addition 
to  plane  flames.   An  important  difference  between  plane  and  non-planar 
flames  should  be  mentioned  here.   For  non-planar  flames  the  spatial 
variable  (radius)  appears  explicitly  in  the  gradient  termiS  v/hile  for 
plane  flames  the  spatial  variables  appear  only  in  the  derivatives. 
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Thus  for  the  simple  cases,  the  space  variable  can  be  eliminated  and 
the  governing  equations  can  be  reduced  to  a  single  equation  \<i±th. 
temperature  as  the  independent  variable  and  the  mass  flux  as  the 
dependent  variable.   An  added  advantage  of  this  is  that  the  range 
of  integration  is  now  finite  in  the  independent  variable. 

Tl-ie  method  of  lines  was  used  to  solve  the  governing  equations 
in  Bledjian's  formulation  [55].   This  consists  of  converting  the  set 
of  partial  differential  equations  to  a  set  of  ordinary  differential 
equations  by  discretizing  all  but  one  of  the  independent  variables 
at  a  time,  the  resulting  equations  being  then  integrated  using 
Runge-Kutta's  fourth  order  method.   Difficulties  were  found  when  the 
set  of  equations  was  stiff.   The  step  size  for  the  next  integration 
cycle  was  estimated  based  on  the  local  truncation  error,  as  calculated 
by  using  Richardson's  extrapolation  formulas.   Calculations  were  carried 
out  for  two  flames:   the  hydrazine  decomposition  flame  and  the  ozone 
decomposition  flame.   The  results  for  the  first  case  agreed  well  with 
Spalding's  calculations  [45]  and  for  the  second  case  good  agreem.ent 
was  obtained  with  experimental  results.   Estimations  of  the  minimum 
ignition  energy  for  the  t^jo  cases  have  also  been  obtained.   For  the 
sake  of  computational  stability,  the  time  step  size  was  restricted  by 

the  condition  (A/pC  )At/(A'^)^  ^   1/30  instead  of  the  standard  stability 

2 
condition  (A/pC  )At(ai|j)   <  1/2  for  the  heat  conduction  equatxon. 
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CHAPTER  IV 
I^ATHEt-IATICAL  ANALYSIS  OF  THE  COMBUSTION  PROCESS 


4.1   Formulation  of  the  Proble 


m 


In  this  chapter  the  case  of  laminar  premixed  flame  propagation 
in  a  stoichiometric  mixture  of  methane  and  air,  initially  at  a  pressure 
of  10  atmospheres  and  at  a  temperature  of  628  degrees  K,  which  approx- 
imately is  the  condition  in  an  Otto  cycle  engine  at  the  time  of  ignition, 
has  been  considered.   The  time-dependent  approach  of  Spalding  [1]  has 
been  used  here.   The  reason  methane  has  been  considered  is  because  of 
the  availability  of  the  reaction  kinetics  data.   The  method  is  general, 
and  if  the  reaction  miechanism  and  kinetic  data  for  the  oxidation  of  any 
other  fuel  are  known,  the  same  method  can  be  used  to  obtain  the  detailed 
flame  structure.   Tiie  variations  of  the  transport  and  thermodynamic 
properties  with  temperature,  pressure  and  chemical  comipositlon  have 
been  accounted  for  as  far  as  the  available  data  permit. 

4.1.1   The  conservation  equations 

The  flame  zone  is  generally  a  few  hundred  mean  free  paths 
thick,  and  so  the  concepts  of  continuum  analysis  may  be  used  here. 

The  general  conservation  equations  for  m:ulti-com.ponent  reacting 
ideal  gas  mixtures  are  given  in  various  references  [2,3],   These 
equations  are; 

(1)     Overall  continuity: 

^P  +  V- (pv)  =  0  (1) 


at 
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(2)     Moment 


um: 


j^  +    (vV)v  = 


(V-p)/p  +  I      Y.f. 
i=l 
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(3)     Energy: 
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(4)     Continuity  of  species: 
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In  addition  to  the  above  conservations  equations,  other  equations 
describing  the  fluxes  of  mass  and  energy,  the  chemical  kinetics  of  the 
system  and  some  constitutive  relations  are  necessary.   Tliese  are  (again 
in  their  general  form) : 

(1)     Diffusion  in  multi-component  mixtures  is  governed  by: 
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and 


i=l 


Y.V. 
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The  derivation  of  equation  (5)  assumes  that  the  medium  under 
consideration  can  be  treated  as  a  dilute  gas.   However,  there  is 


39 


experimental  evidence  that  this  result  also  applies  quite  well  in 
dense  gases  and  liquids  [3,  p.  718]. 


(2)     Heat  flux  vector: 
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The  same  comment  as  above  applies  to  this  equation  also. 

(3)     Chemical  kinetics  for  a  set  of  M  elementary  reactions  between 
N  species: 
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(4)     Thermally  perfect  gas: 


P  =  pR  T/W 
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(9) 


w  =  y  x.w 


1  =  1 


(5)     Caloric  equation  of  state; 
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Implicit  in  these  equations  is  the  assumption  that  every  point 
in  the  flow  field  has  a  thermodynamically  definable  state,  i.e.,  it  is 
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possible  to  attach  significance  to  the  concepts  of  local  temperature, 
pressure,  density  and  chemical  composition  (see  Sec.  4.4  for  un- 
certainties in  the  definition  of  temperature). 

4.1.2   Major  assumptions  for  analysis 

The  general  equations  are  now  simplified  for  the  case  of  one- 
dimensional  laminar  flame  propagation  in  spherical  coordinates  using 
the  following  assumptions. 

(1)  The  kinetic  energy  associated  with  the  motion  of  the  gases  is 
small  compared  to  the  heat  of  the  reaction  and  may  be  neglected.   Also, 
since  the  velocities  and  velocity  gradients  are  small,  viscous  effects 
may  be  neglected.   Thus  the  momentum  equation  can  be  replaced  by  the 
assumption  that  the  pressure  is  constant  through  the  flame.   In  fact, 
solving  the  momentum  equation  for  the  pressure  change  across  the 
flame  yields: 

r   V  ^  s 

^P  =  p.  -  P   =  p  Y   1  -  TT  M   ;    M  -  --  (12) 

'^b    u    u       Vu      u   a 
*-     u'^  u 

From  experimental  observations  it  is  seen  that  in  general 

M  <  0.02,  V,  /V   <  15,  and  Y  <   1.667,  which  yields  Ap  ^  0.01  p  . 
u  b   u  u 

Hence  for  ordinary  flames  the  process  may  be  considered  essentially 
isobaric.   This  assumption  v7ould  not  hold  for  detonation  waves. 

(2)  Radiative  heat  transfer  is  negligible.   Although  relatively 
high  temperatures  (~  2500°K)  are  encountered  in  flames,  radiation  is 
not  important  since  the  gases  have  low  emissivity  and  absorptivity  and 

so  the  radiative  heat  transfer  during  the  reaction  time  is  small  compared 
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to  the  total  energy  flux.   With  flames  containing  solid  carbon 
this  is  generally  not  the  case. 

(3)  Ihe  Soret  and  Dufour  effects  are  negligible. 

(4)  Any  diffusion  caused  by  pressure  gradients  is  negligible.   Since 
assumption  (1)  already  considers  the  process  to  be  Isobaric,  this 
particular  assumption  is  not  necessary.   This  term  may  be  significant 
in  strong  detonation  waves  wherein  tremendous  pressure  gradients 

may  be  established. 

(5)  External  body  forces  are  neglected.   For  very  hot  flam.es  wherein 
ionization  occurs  to  a  significant  degree,  this  term  may  have  to  be 
considered  if  an  external  electric  field  is  also  present;  the  external 
force  on  an  ion  is  equal  to  the  product  of  the  ionic  charge  and  the 
local  electric  field  strength.   This  gives  rise  to  what  is  termed 

as  'forced  diffusion'. 

(6)  The  chemical  kinetics  is  assumed  to  follow  Arrhenius  equation. 
Thus  : 


kf   =  A^T^  exp(-r^/RJ)  (13) 


Tlie  backward   reaction   rate    constant    ku      is    then   obtained   from 

■^k 

a  knowledge  of  kp   and  the  equilibrium  constant  IC-^  : 

^fk 


There  is  considerable  controversy  [4,5]  regarding  the  validity 
of  equation  (14);  it  will  nevertheless  be  used  in  the  absence  of  a 
more  accurate  relation. 
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(7)   The  gas  mixture  is  thermally  perfect,  i.e.,  it  obeys  the  ideal 
gas  equation.   Nitrogen  is  a  major  constituent  (>  70%)  of  the 
mixture  and  it  closely  approximates  thermally  perfect  behavior  for 
the  conditions  considered  in  this  study.   Other  constituents  are 
present  in  relatively  smaller  concentrations,  so  as  not  to  affect 
the  overall  thermally  perfect  behavior. 

The  equations  finally  simplify  to: 


(1)     Overall  continuity: 


■^  +  -5-  TT-  (r  pv  )  =  0 
dt     2  dr      r 
r 


(15) 


(subscript  r  denotes  the  radial  direction) 


(2)     Energy  equation: 


9T  ,  dl 
"57"  "f"  V  3— 
at  r  dr 


119 

pC   2  3r 
P  r 


rh^l    -^  I      Y.V.   -^ 
dr     C   .^,   iirdr 
^  p  1=1 


pC   /    IX 
p  1=1 


(16) 


(3)      Species  continuity: 


3Y.       3y.    w. 

i+v  ^=--i-i.i-^  (r^pY.V.  )  ,   i=l,2, 
dt     r  or    p    p    2  dr      1  ir 


N   (17) 


(4)     Diffusion  equation; 


3X.     N   /-X.X  - 


^-      I 


j=l 


^j  ' 


(V.   -  V.  )  ,    i=l,2 
jr    ir   '       ' 


.N 


(18) 
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N 
and  y  Y.V.   =0  (iQ) 

J=l    -^ 


(18)  represents  a  set  of  (N-1)  linearly  independent  equations,  since 

N  dX. 

.  ^,  dr 
1=1 

The  unknowns  in  these  equations  are:   Y. ,  V.  ,  T,  v  .   X.  and 

X   ir      r    1 

p  are  determined  once  Y.  and  T  are  kno\^m  respectively,  using  the 
following  two  equations: 


Y.  /W. 
X.  =   i^  ^  "- ,    i=l,2, N  (20) 

I      (Y,/W  ) 
j  =  l   ^   ^ 

p "  W 
and  p  =  g-Y   (the  ideal  gas  equation)  (21) 

u 


4.1.3   Boundary  conditions 


The  next  step  is  to  set  up  the  boundary  conditions  at  the  hot 
and  cold  ends  tor  the  set  of  partial  differential  equations  given 
above. 

a.   Ihe  cold  boundary  difficulty.   This  has  been  discussed 
considerably  in  the  literature.   Basically  the  difficulty  arises  from 
the  incompatibility  of  the  energy  conservation  equation  with  the 
boundary  condition  it  is  subjected  to  at  the  cold  end.   The  model  for 
the  flame  which  assumes  both  'steady  state'  and  'adiabatic'  conditions 
does  not  represent  truly  the  physical  processes  occurring  and  these 
conditions  contradict  each  other  for  a  realistic  reaction  rate  law. 
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wherein  the  reaction  rate  at  the  cold  boundary  is  small  but  finite 
[6].   To  satisfy  the  conditions  of  both  adiabatic  and  steady  state 
requires  any  one  of  several  assumptions.   The  assumptions  are 
mathematically  equivalent  but  have  different  physical  meanings. 
Also  the  validity  of  these  assumptions  on  physical  grounds  is  open 
to  question,  although  these  assumptions  have  been  used  widely  to 
solve  the  flame  equations.   The  assumptions  are: 

(1)  This  relaxes  the  adiabatic  condition  at  the  cold  boundary 
by  alloxving  the  existence  of  a  small  but  finite  heat  sink  at  the 
boundary  in  the  form  of  a  flame  holder.   It  was  introduced  originally 
by  Hirschfelder  et  al.  [3].   This  was  done  because  a  solution  does 
not  exist  when  the  very  small  heat  loss  at  the  cold  boundary  is 
neglected;  the  neglect  of  heat  transfer  (due  to  lateral  conduction 
and  radiation)  at  other  points  in  the  flame  zone  does  not,  however, 
affect  the  existence  of  a  solution.   Mathematically,  this  is  equivalent 
to  relaxing  the  condition  (8T/3r)  =  0  at  the  cold  boundary,  and  instead 
allowing  for  (3T/3r)  =  e,  where  e  is  a  small  positive  quantity.   The 
value  of  c  determines  the  'strength'  of  the  flameholder  in  terms  of 

the  heat  abstracted  by  it  and  has  to  be  greater  than  a  certain  critical 
value  to  permit  the  existence  of  a  solution.   This  assumption  is  physically 
unrealistic  for  a  freely  propagating  flame  which  does  not  have  a  flame- 
holder  or  heat  sink  propagating  xjith  it. 

(2)  The  second  assumption  compromises  another  physical  aspect 
by  allowing  for  the  existence  of  a  temperature  (greater  than  the  initial 
temperature  of  the  reactants)  below  I'/hich  the  heat  release  due  to 
chemical  reactions  is  identically  zero.   This  tem/perature  has  widely 
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been  referred  to  as  the  'ignition  temperature'  [7].   The  choice  of 

this  temperature  is  again  not  arbitrary  but  depends  on  the  particular 

flame  being  studied  and  it  can  be  determined  from  the  governing 

equations.   It  permits  the  adiabatic  assumption  (9T/9r)  =  0  at 

T  =  T   to  be  satisfied:  however,  now  this  occurs  at  r  =  °°,  as 
u  '  ' 

required  for  a  strict  mathematical  model  [6]. 

The  use  of  the  ignition  temperature  concept  can  be  explained 
by  relaxing  the  adiabatic  condition.   It  is  assumed  that  the  heat 
loss  from  the  cold  boundary  and  the  heat  released  by  chemical 
reactions  (not  identically  zero)  at  the  cold  boundary  reach  an 
equilibrium  at  the  temperature  T  .   Further,  assuming  the  rate  of 
the  chemical  reaction  to  be  extremely  slow  (the  initial  mixture  is 
in  a  metastable  state)  ,  the  change  in  concentration  v;ith  time  at  the 
cold  boundary  can  be  neglected.   Thus  it  is  as  though  no  reaction  has 
taken  place,  the  heat  released  at  the  cold  boundary  having  been 
conducted  away. 

(3)   ITie  third  assumption  involves  the  use  of  truncated 
kinetics — v/herein  the  reaction  rate  goes  to  zero  as  the  temperature 
approaches  its  cold  boundary  value.   Tliis  has  been  shown  to  be 
equivalent  to  the  use  of  an  ignition  temperature  which  is  very  close 
to,  but  not  equal  to  the  cold  boundary  temperature.   Again,  this 
assumption  is  not  physically  possible.   However,  for  most  systems 
of  interest,  the  reaction  rate  at  the  cold  boundary  is  so  small  as 
to  be  immeasurable,  and  so  it  can  be  considered  to  be  essentially 
zero.   'ibe  reaction  rate  law  is  then  written  as 
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k^  =  at"  exp(-E/RT^T)                        (22) 

T  -  T 

where               x  =  "                                (23) 

b  u 


as  compared  to  the  actual  law. 


k^  =  at""  exp(-E/RT)  (24) 


The  factor  T  varies  from  zero  at  the  cold  boundary  to  unity  at 
the  hot  boundary.   The  change  in  magnitude  of  the  reaction  rate  due  to 
this  modification  becomes  smaller  as  the  temperature  approaches  its 
hot  boundary  value  and  this  is  where  most  of  the  chem.ical  reaction 
and  its  associated  heat  release  occur. 

In  the  calculations  carried  out  here  for  the  methane-air 
flam.e  with  an  initial  temperature  of  628  degrees  K,  it  was  found 
that  the  chemical  reaction  rate  at  this  temperature  was  lower  than 
the  sm.allest  number  that  can  be  handled  by  the  IBM  370/165  computer 
[i.e.,  exp(-176)].   In  fact,  the  reaction  rates  were  ignored  at 
temperatures  belox\'  650  degrees  K  for  all  the  elem.entary  reaction 
steps,  and  above  650  degrees  K,  the  reactions  for  which  E/RT  T  was 
greater  than  120.0  were  ignored.   The  condition  E/RT  T  >  120.0  was 
used  instead  of  E/RT,  T  >  176.0  to  avoid  the  occurrence  of  a  number 

D 

smaller  than  exp(-176)  later  during  the  com.putations . 

Again,  when  integrating  from  the  hot  to  the  cold  boundary, 
the  integration  was  stopped  arbitrarily  for  I  (F^"  '^ic-i  ^/"^ic  I  ~  '-'•'^O-^; 
(F  =  dependent  variable;  temperature  and  mass  fractions)  the  subscript 
'k'  refers  to  the  grid  point  in  the  r-direction.   This  had  to  be  done 
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since,  when  carrying  out  numerical  calculations,  the  condition 
(9F/3r)  =  0  can  strictly  be  achieved  only  at  r  =  =°. 

b.   Hot  boundary.   Iftien  only  the  initial  (unreacted)  and  final 
(reacted)  states  of  the  gas  are  considered,  it  is  seen  that  the  gas 
has  passed  through  the  flame  zone  without  any  net  gain  or  loss  of 
heat.   All  that  has  occurred  is  a  change  in  composition  due  to 
chemical  reaction  and  all  the  heat  of  the  reaction  has  gone  to  raise 
the  temperature  of  the  gas.   This  adiabatic  condition  is  not  fulfilled 
for  f]ames  which  lose  significant  amounts  of  heat  by  other  processes, 
e.g.,  radiation,  heat  losses  to  f  lameholders ,  etc.   Hox^rever,  within 
the  fram.ework  of  the  assumptions  used  in  this  analysis,  the  adiabatic 
assumption  is  satisfied.   Again,  the  adiabatic  condition  is  not 
satisfied  for  any  elem.ent  of  gas  within  the  flame  zone,  since  it 
loses  or  gains  heat  by  conduction  and  diffusion  of  species.   ITierefore, 
for  a  constant  pressure  flame,  since  the  heat  transfer  equals  the 
change  in  enthalpy,  the  enthalpies  (per  unit  mass  of  mixture)  of  only 
the  initial  and  the  final  states  are  equal. 

For  the  special  case  of  unity  Lexvis  number,  the  constant 
enthalpy  condition  is  satisfied  throughout  the  flame  zone.  Since 
the  heat  fluxes  due  to  diffusion  and  thermal  conduction  are  equal 
in  such  a  case,  any  element  of  gas  within  the  flame  zone  does  not 
suffer  any  net  gain  or  loss  of  heat  and  so  goes  through  the  flam.e 
zone  at  constant  enthalpy. 

Tne  condition  of  equal  enthalpy  at  the  hot  and  cold  ends  is 
utilized  to  obtain  the  boundary  condition  at  the  hot  end.   Since  the 
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gas  mixture  is  in  chemical  equilibrium  at  the  hot  end,  the  concentra- 
tions of  the  various  species  are  obtained  by  simply  calculating  the 
equilibrium  concentrations  corresponding  to  an  isenthalpic  combustion 
process,  starting  with  a  mixture  of  given  concentrations  and  state 
properties  at  the  cold  end,  i.e., 


h(T^,p)  =  h(T^,p)  (25) 


which  is  known  from  the  initial  conditions  at  the  cold  end.   Thus 


N 
h(T  ,p)  =  I    (h  Xi  )  (26) 

U  .   -,    J-    11  J-— J- 


i=l 


where  the  h.  are  obtained  from  the  JANAF  Tables. 

X 

Considering  a  Taylor  series  expansion  about  Tjj  ,  the  kth 
estimate  for  the  value  of  the  burnt  gas  temperature,  and  neglecting 
higher  order  terms. 


h(T^,p)  =  h(T^  ,p)  +  (T^^  -  T^^) 


3h 


3TJ 
'   P 


(27) 


where  the  partial  derivative  is  evaluated  at  (T^^  ,p)  and 

k 


h(Tb   p)  =   I   Vi.^T^T.  ^2^) 


h(T,  ,p)  -  h(Tb.  ,p) 


"".  =  — MV) —  ''" 


and  Th,  ,,  .  Th.  +  ATu,.  (30) 


Tb^+1  -  Tb^  +  AT^^^ 
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A  first  estimation  of  the  temperature  at  the  hot  end  is 
obtained  from: 


T  =  T  +     fuel 
b    u    C  (T  ) 
P   u 


(31) 


where   Ah  =  heat  of  combustion  (kcal/gm.mole) 


C  (T  ) 
P   u 


specific  heat  at  constant  pressure  of  the  mixture  at 
temperature  T   (kcal/gm.mole  °K) 


1=1  ^ 


(32) 


The  C    are  again  obtained  from  the  JANAF  Tables. 
^i 

Tlie  equilibrium  concentrations  of  the  species  are  calculated 
at  the  estimated  value  of  the  burnt  gas  temperature  by  the  method 
given  in  Appendix  I.   Then  equation  (30)  is  used  to  obtain  the  next 
estimate  of  this  temperature  and  the  process  is  repeated  till 
convergence  is  obtained  to  within  a  predetermined  accuracy: 


^Tbv 


E  =   0.001  was    used   in   the    calculations. 

This    determines    the  boundary    conditions    at    the  hot   end: 


(33) 


T   =    T, 


X.    =   X. 
1  ^b 


i=l,2 


(34) 
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4.2   Method  of  Solution 

An  explicit  finite  difference  method  was  used  for  numerically 
solving  the  partial  differential  equations  (15-19),  and  obtain  the 
profiles  of  the  state  properties  and  the  species  concentrations 
through  the  flame  zone.   A  scheme  involving  a  symmetrically  staggered 
mesh  in  the  r-t  plane  was  used.   Tliis  method  has  been  used  by  Nash  [8] 
for  the  solution  of  the  three-dimensional  boundary  layer  equations. 

The  staggered  mesh  scheme  uses  a  set  of  secondary  points  in 
addition  to  the  primary  collocation  points,  the  former  lying  at  the 
centers  of  the  rectangles  formed  by  the  latter,  as  shown  in  Fig.  4.1. 

In  this  method,  values  of  the  dependent  variables  are  at  first 
calculated  at  the  secondary  points  [at  time  (t+At/2)]  using  the  values 
associated  with  the  primary  points  at  time   t  .   Then  using  values 
associated  with  both  tVie  primary  points  at  tim.e  t  and  the  secondary 
points  at  time  (t+At/2),  values  of  the  dependent  variables  at  the 
primary  points  at  time  (t+At)  can  be  obtained.   This  is  more  accurate 
than  proceeding  directly  from  time  level  t  to  level  (t+At). 

The  difference  equations  for  the  two  steps  are  readily  obtained 
by  double  Taylor  series  expansions  about  the  appropriate  points,  and 
are  [8]: 


f 


^c   =  2  ^^'a  ■"   ^> 


AAi 


'3f" 
3r 

A 

BJ 

■  +  -^  At   i 
4 

+ 

A 

[atj 

(35) 
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hH 


Ar/2  —1  h-  t 


Fig.  4.1(a).   Staggered  mesh  scheme.   r-t  plane 


t  . 


O  primary  points 
n  secondary  points 


E(i,k+1) 


Fig.  4.1(b).   Enlarged  view  of  staggered  mesh 


lW»»««WW"l»WrW>M!IBW<i^^W«a*|iWpiAr||ipqilMW»<IM>i>waL!»W»MWi^ 
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Fj,  =  UF^  +  2  (l-y)(Fj3  +  F^)  +  X(l-y)Ar 


(d7 

D 

[arj 

+  ^  (i+y)At 


+  e. 


(36) 


Here  F  denotes  any  one  of  the  dependent  variables.   The  quantities 
(3F/8t)  are  obtained  from  the  governing  partial  differential  equations 
(16)  and  (17),  the  r-derivatlves  in  these  equations  and  in  equations  (35) 
and  (36)  being  expressed  in  finite  difference  form.   A  fourth  order  dif- 
ference formula  is  used  for  this  purpose: 


3F 


dr 


(F 


)  - 


(F_. 


-  F. 


i,k 


12Ar  '  i-2,k    i+2,k'    3Ar  '  i-l,k    i+1 ,k 


+  30  ^^^-^  TT 


) 


(37) 


In  equations  (35)  and  (36),  e   and  e„  denote  the  discretization 
or  truncation  errors,  and  contain  higher  order  terms  of  the  double 
Taylor  series  expansions. 


A         T>%/A  ^2  9^F    1  ,.  ,2  3^F  _^  X  ,  ,,  -2   9^F  ^ 
e,  =  (s-A)(Ar)  _  -  -  (At)      +-At(Ar)  -~—^  + 

or  dt  dtdr 


(38) 


r-,       ^A         ^s,A    ^2  9^F    1  .    ...  .2  3^F 
e^  ==  (l-y)(g  -  A)(Ar)  — ^  "  8  (l-l^)('^t)  —j 

dr  3t 


3 

+  i  (4A-4Ay  +  M)At(Ar)^  -^-^  + 


(39) 


A  and  y  are  constants  which  govern  the  accuracy  and  the 
stability  of  the  computation.   The  pair  of  values  A  =  1/8  and  y  =  0 
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serve  to  minimize  e„  and  also  assure  stability  [8],   With  these  values 
the  two-step  scheme  reduces  to  a  one-step  scheme: 


^C  =  I  (^A  +  ^>  "-   I  ^^ 


3F 


8F 


4 


f9F^ 

+ 
A 

[9tJ 

B. 

(40) 


^E  =  i(^D  +  ^c)+l^- 


D 

3F 
l3rj 

4.. 


f3F^ 

+ 
D 

[BtJ 

C 

(41) 


4.2.1   Initialization 

S-shaped  profiles  for  the  dependent  variables  Y.  and  T  were 
assumed  to  start  the  integration.   Second  order  curves  were  chosen  and 
these  were  of  the  form: 


F(r) 
F(r) 

F(r) 
F(r) 


=  F 


r   d   r 


F  +  (F,  -  F  ) 
u     b    u 


F   -  (F  -  F  ) 
b      b    u 


(r  -  r)' 


26^ 


9 

26^ 


\    '  ''b  ^  ^ 


r   >  r  >  r,  +  6 
u        b 


r,  +  6  >  r  >  r, 
b  b 


(42) 


An  arbitrary  value  of  6  =  —  (r  -  r,  )  =  0.005  cm.  was  chosen. 
It  is  immaterial  what  value  of  6  is  chosen  since  the  final  profiles  are 
independent  of  the  starting  ones. 

4.2.2  Details  of  computation 

The  initial  profiles  having  been  fixed,  the  integration  process 


can  now  be  started.   Expressing  (9X./8r)  in  finite  difference  form  in 
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accordance  with  equation  (37),  the  diffusion  equations  (18)  and  (19) 

reduce  to  a  set  of  N  simultaneous  algebraic  equations  with  the  N 

diffusion  velocities  V.   (j=l ,2 , . . . .N)  as  the  unknowns.   These 

Jr 

equations  were  solved  using  Gauss-Jordan  elimination  with  maximum 
pivot  strategy.   Thus  the  diffusion  velocities  were  obtained  at  each 
grid  point  for  a  given  time  level  t. 

Before  proceeding  further,  tx^ro  other  methods  for  the  determina- 
tion of  diffusion  velocities  will  be  mentioned.   These  methods  are 
approximate,  but  involve  lesser  computation  than  the  method  mentioned 
above,  and  may  therefore  be  used  where  conditions  warrant  their  use. 

(1)   In  flames  involving  air  as  oxidant,  one  component  (namely,  nitrogen) 

is  present  in  great  excess  so  that  it  is  possible  to  simplify  the  problem 

by  regarding  each  of  the  other  species,  one  at  a  time,  as  forming  a 

binary  mixture  with  the  main  component.   In  general,  the  effective 

binary  diffusion  coefficient  D.   for  the  diffusion  of  species  i  in  a 

im  ^ 

mixture   m  is    given  by    [9]: 


y       (1/nD. .)(X.N.    -   X.N.' 
1  1=1  iJ         J    1  ^    J 


(^3) 
nD.  N 

1™  M  ^r  V  TIT 

N.    -  X.       >       N. 


For  the  diffusion  of  trace  component  i=l,2,3,...N  (ir j ) ,  in  nearly  pure 
species  j , 


D.   ~  D. . 
im    ij 


and  so  the  m.olar  fluxes  and  hence  the  diffusion  velocities  arc  obtained 
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by   solving   the  set   of   algebraic  equations: 


dX.  N 

N.    =  -  nD.      —  +  X.  y      N.                                                 (44) 

1                 xm   dr             1  .    ,      i                                                 ^      ' 

N.    =  n.Vi^  (45) 


(2)   This  approximation  may  be  used  when  species  i  is  still  a  trace, 
but  there  is  no  single  component  which  may  be  considered  to  be  in 
excess  compared  to  all  the  others.   For  such  a  case  the  diffusion 
velocity  of  the  trace  species  is  obtained  from: 


V   = -L _3.  (46) 

X.   )   X./D. 


dr 


v/hich  shows    that    the   quantity 


-1 


is  an  effective  diffusion 


coefficient  for  the  species  i  and  the  mixture.   The  above  equation  is 
not  valid  for  the  diffusion  of  a  species  present  in  significant  amounts. 
Such  species  may  be  treated  in  two  ways.   Often  a  species  which  is 
present  in  excess  is  not  significantly  affected  by  the  chemical  reaction 
(e.g.,  nitrogen)  so  that  its  concentration  is  m.ore  or  less  constant 
through  the  flam.e  zone  and  hence  its  diffusion  velocity  may  be  neglected. 
On  the  other  hand,  for  species  which  are  chemically  reactive,  and  not 
present  as  a  trace,  the  diffusion  velocity  may  be  computed  by  difference 
from  equation  (19),  which  is 

N 

I   Y.V.   =  0 
.^,   X  ^r 
x=l 
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Getting  back  to  the  main  problem,  it  is  seen  that  the  mass- 
average  velocity  v  is  as  yet  unknown,  and  a  suitable  profile  for  it 
was  assumed  subject  to  Vj-  =  0,  i.e.,  the  unburnt  gases  are  at  rest. 
Tliis  need  be  done  only  for  integration  over  the  first  time  step;  for 
succeeding  time  steps,  the  values  of  v  from,  the  preceding  time  step 
were  used  for  a  first  approximation.   With  this  assumed  profile  for 

v^,  (9T/9t)  and  (8Y./3t)  (i=l,2, N)  were  calculated  at  each  grid 

point  for  a  given  time  level  t,  using  equations  (16)  and  (17).   The 
derivatives  in  r  occurring  in  these  equations  Xizere  approximated,  again 
by  the  fourth  order  difference  form.ula  (37);  values  of  Vj   (j=l  ,2 , .  .  ,  .N) 
as  calculated  above  were  used,  and  the  JANAF  Tables  were  used  to  obtain 
the  values  of  the  state  properties  h  and  C  .   The  V7.  (i=l  ,2 , .  .  .  .  N)  were 
calculated  according  to  equation  (8)  and  modified  by  truncating  the 
reaction  kinetics  as  discussed  in  connection  with  the  cold  boundary 
difficulty. 

Using  the  ideal  gas  equation,  it  may  be  shown  that: 


*■  1=1  1 


and  thus  the  overall  continuity  equation  may  be  rev/ritcen  in  the  form: 


1 
V   = 


r     2 

Pr 


""   r   N     3Y.    ^  ,„, 

'   ^  1=1  1 


r^dr  (48) 


r 

u 


(for  given  time  t) 
satisfying  v^.  ^  =  0. 


■u 
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The  values  of  (3T/3t)  and  (9Y./8t)  calculated  above  with  an 
assumed  profile  for  v  were  now  used  in  equation  (48)  to  obtain  better 
values  for  v  ,  which  were  then  used  to  obtain  improved  values  of 
(8T/9t)  and  (9Y./3t).   The  iteration  can  be  terminated  when  the 
v  values  converge  to  within  a  predetermined  accuracy. 

Actually  the  values  of  (9T/8t)  betv;een  two  successive  iterations 

were  compared  for  checking  convergence.   This  was  done  since  (9T/9t)  is 

used  to  determine  the  temperature  at  the  next  time  level  (equations  40, 

41)  and  so  convergence  of  (3T/3t)  is  of  greater  importance  than  that 

of  V  which  is  not  used  any^'/here  else  in  the  calculations.   The  con- 
r 

vergence  of  (3Y./3t)  (i=l ,2 , . . . .N)  may  also  be  tested,  but  these  were 
not  found  to  vary  significantly  and  the  convergence  of  (3T/3t)  generally 
assures  the  convergence  of  (3Y./3t)  for  all  i.   A  convergence  criterion 
of  the  type 


At 


3T 
3t 


3T 
3t 


k-lj 


^  e 


(49) 


v;as  used  with  £  =  0.0005.   This  assures  the  errors  in  the  calculated 
temperatures  due  to  uncertainty  in  the  term  (3T/3t)  to  be  less  than 
e/2  each  time  either  equation  (40)  or  (41)  is  used.   The  subscripts 
in  the  above  convergence  condition  refer  to  the  iteration  number  and 
this  condition  is  to  be  satisfied  at  each  grid  point  in  the  flam.e 
zone  before  proceeding  further. 

Equation  (48)  involves  quadrature;  this  ^^7as  done  using  the 
trapezoidal  rule  with  the  same  space  interval  Ar  as  used  for  the 


other  computations.   Simpson's  rule  along  with  Newton's  3/8  rule  may 
also  be  used  (IKl  subroutine  QSF),   The  difference  between  the  two 
was  found  to  be  negligible  since  the  space  interval  Ar  is  small. 
Hence  the  trapezoidal  rule  which  needs  less    computer  running  time 
was  used. 

The  term  involving  (3Y./9t)  in  equation  (48)  has  only  a  small 
influence  on  the  final  profile  of  v  and  may  be  neglected  to  a  first 
approximation.   This  would  amount  to  making  the  assumption  that  the 
average  molecular  weight  of  the  gas  mixture  is  constant  through  the 
flame  zone.   This  is  approximately  true  for  methane-air  flames  in 
which  a  single  constituent,  nitrogen,  accounts  for  more  than  70%  of 
the  concentration  at  any  point  in  the  flame  zone.   For  the  sake  of 
completeness,  this  term  was  included  in  the  calculations  carried  out 
here,  although  its  effect  is  small. 

The  iterative  scheme  used  here  for  the  solution  of  equation  (48) 
is  similar  to  Picard's  method  of  successive  approximations  for  the 
solution  of  a  Volterra  integral  equation  [10].   This  is  seen  from  the 
analysis  given  below. 

Equations  (16)  and  (17)  may  be  written  as; 


■^=  -  V  |^+G(r)  (50) 

dt  r    or 


dY.  3Y. 

■^='\ir-'\^'^  '    ^==^'2, N  (51) 


where 
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G(r)  =  - 


1   1  d 


pC   2  dr 
P  r 


2,  8T 

r  A  -=r— 

dr 


N        8h.     T    N 
TT-     I      Y.Vi  -^  -  -^  I     h.w.   (52) 
p  1=1  p  1=1 


and 


«i(^>=p--^\l7(^Wir)  '    -^'2' ^ 


(53) 


At  a  given  instant  of  time  the  right  hand  sides  of  (52)  and  (53) 
are  functions  of  r  only. 

Equation  (48)  now  reduces  to: 
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Again,  for  a  given  instant  of  time 
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(54) 


N   ^   3Y.    ^  .^ 
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N  H.  (r)    „,  . 

1=1     1 


pr   =  A^(r) 


(55) 


where  A   A   and  A  are  functions  of  r  only.   So  (54)  finally  reduces 
to: 


v^(r) 


A3(r)  . 


A^(r)  [v^(r)A^(r)  +A2(r)]dr 


(56) 
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This  is  nothing  but  a  Volterra  integral  equation  of  the  second 
kind,  and  the  iterative  technique  described  above  is  Picard's  method  of 
successive  approximations.   The  first  approximation  to  v  (r)  in  this 


method  is  given  by: 


^^"^  =  a7(7)  J   -3-^-2 


A^(r)A^(r)dr 
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Having  finalized  the  values  of  the  time  derivatives  (3T/3t)  and 
(oY./3t)  at  a  given  time  instant  t,  equation  (40)  was  used  to  proceed 
to  the  secondary  points  at  the  time  level  (t+At/2)  and  by  a  similar 
process  the  values  of  the  dependent  variables  at  time  level  (t+At)  were 
obtainad  using  equation  (41).   Ihe  mole  fractions  were   then  calculated 
from  the  mass  fractions  using  equation  (20)  and  the  m.ean  gas  density 
was  obtained  from  a  knowledge  of  the  temperature  and  the  m.ean  molecular 
weight  using  the  ideal  gas  equation.   Calculations  in  the  r-direction 
at  the  cold  boundary  were  stopped  at  any  given  time  level  when  it  was 
found  that  the  relative  changes  in  the  temperature  and  the  species 
mass  fractions  were  less  than  a  predetermined  value,  i.e.. 
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F  -  F 
k    k-1 


\ 


<   e  (58) 


e  =  0.005  was  used  in  the  calculations.   This  process  X'/as  repeated 
until  the  dependent  variables  converged  to  their  steady  state  values. 
The  transport  properties  in  the  conservation  equations  were  determined 
(Appendix  II)  at  each  grid  point  every  75  time  steps;  the  error 
introduced  by  not  calculating  the  new  values  of  the  transport 
properties  at  every  time  step  would  cause  small  errors  in  the  transient 
state  solutions,  but  would  not  affect  the  steady  state  solution.   However, 
since  the  estimation  of  the  transport  coefficients  is  quite  often  approx- 
imate, the  error  introduced  by  not  calculating  these  coefficients  at 
every  time  step,  might  in  fact  be  lower  than  the  error  involved  in  the 
calculations  based  on  the  approximate  theories. 

The  last  step  involves  the  computation  of  the  laminar  flam.e 
velocity  based  on  the  unbumt  gas  density.   This  was  done  by  considering 
a  control  volum.e  which  is  large  compared  to  the  thickness  of  the  flame 
zone,  and  equating  the  net  rate  of  production  of  species  i  due  to 
chemical  reactions  to  the  rate  of  outflow  of  species  i.   Tlie  result  is: 


^b 

w.dr 

1 

r 

p  S  =  ^ ^T—  (59) 

u  u   Y.  -  Y., 
lu   lb 

The  values  of  the  flame  velocity,  S  ,  as  calculated  from  the 

J  J   u 

above  equation  converged  as  the  state  properties  and  species  concentra- 
tions approached  their  steady  state  values  only  when  species  i  was 
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considered  to  be  one  of  the  major  species  (i.e.,  CO^ ,  H2O,  CH^ ,  "^2^' 

however,  for  each  one  of  the  species  a  different  value  of  the  flame 

velocity  was  obtained,   Ihe  reason  for  this  are  the  unreliable  kinetic 

data.   Further,  this  type  of  behavior  has  been  observed  before  [11]. 

The  flame  velocity  calculations  based  on  the  other  species  did  not 

converge  to  a  steady  value,  presumably  because  the  term  (\\^  -  ^^^) 

involved  in  the  denominator  of  equation  (59)  is  small  for  these 

species  and  sm.all  uncertainties  in  this  term  could  cause  large  errors 

in  the  calculation  of  S  .   Tliese  values  might,  however,  converge  if 

u 

computations  are  carried  on  further.   In  the  calculations  reported 

in  Ref.  12,  it  was  found  that  the  flame  velocity  based  on  the  unstable 

species  needed  more  time  steps  to  converge  to  a  steady  value. 

4.2.3   Step  size  and  stability 

Tne  S-shaped  profiles  used  initially  were  at  least  qualitatively 
correct  for  the  major  species  (CO^ ,  CH^,  0^ ,  and  H2O),   However,  this 
V7as  not  so  for  the  intermediate  species  which  have  a  maximum  concentra- 
tion somewhere  in  the  flame  zone  and  decrease  towards  both  the  hot  and 
the  cold  boundaries.   For  this  reason  initial  time  steps  were  restricted 
to  a  very  small  value  until  the  profiles  of  the  intermediate  species 
approached  their  correct  form,  after  which  larger  step  sizes  could  be 
used.   A  starting  step  size  of  At  =  l.OE-08  sec.  v/as  used  and  this  was 
doubled  every  50  steps  until  it  reached  a  maxim.um  value  of  32.0E-08  sec. 

The  step  size  had  to  be  restricted  to  this  maximum  value  to  satisfy  the 

2 
stability  condition  (X/oC  )At/(Ar)   <  1/2.   Since  A,  p,  and  C  vary 

p  P 

from  point  to  point,  the  combination  of  values  of  A,  p,  and  C  which 
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would  give  the  least  upper  bound  for  At  was  used  to  show  that 

At    ^  1.28(Ar)"  which  gives  At    ^  32.0E-08  sec.  with  Ar  =  5.0E-04  cm. 
max  max 

Thus  this  explicit  scheme  is  conditionally  stable. 

The  use  of  these  values  of  At,  however,  caused  instability  in 
the  integration  of  the  term  involving  the  chemical  rate  of  production 
of  species.   So  this  term  was  integrated  separately  with  an  appropriate 
value  of  step  size.   The  integration  can  be  done  by  Runge-Kutta's  fourth 
order  method  for  a  set  of  first  order  differential  equations.   The 
equations 

w.    dn. 

W-=di^=  fi^^'^2' V'    i=l'2, N       (60) 

1 

?v 
form  a  set  of  what  is  termed  as  'stiff   equations,  and  such  equations 

need  very  small  step  sizes  and  hence  prohibitively  long  computer  running 
times  with  Runge-Kutta's  method. 

Gear  [13]  has  developed  an  algorithm  for  dealing  with  stiff  equa- 
tions, and  this  was  used  in  the  integration  of  the  reaction  rate  equations 
(details  of  the  method  in  Appendix  HI).   Instead  of  starting  the  integra- 
tion by  a  first  order  method,  as  proposed  by  Gear,  a  higher  method  was 
used.   lliis  necessitated  the  calculations  of  up  to  six  derivatives  of  n., 
since  the  maximum  order  for  x-zhich  the  algorithm  is  devised  is  six. 
Generally,  it  was  found  in  this  problem  that  a  starting  fourth  order 
m.ethod  resulted  in  the  maxim.um  value  of  the  step  size.   The  order  for 
subsequent  time  steps  was  decided  by  the  program  itself.   The  following 
equations  were  used  to  obtain  the  derivatives. 


So  called  because  this  behavior  was  first  noticed  in  connection 
with  servomechanisms  having  a  tight  coupling  between  the  driving  and  the 
driven  member. 
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where 
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The  calculation  of  these  derivatives  requires  a  knowledge  of 
the  temperature  to  determine  the  values  of  the  rate  constants  involved 
in  these  equations.   The  temperature  was  assigned  its  value  at  the 
beginning  of  a  time  step  of  the  main  integration  and  regarded  as 
constant  over  this  main  time  step.   The  computations  might  be  repeated, 
however,  using  a  suitable  average  value  once  an  estimate  has  been 
obtained  for  the  temperature  at  the  end  of  the  time  step.   Instead  of 
using  a  simple  arithmetic  average  T  =  y  (T  +  T  ),  an  average  of  the 
nature  exp  (-E/T  R)  =  ~  [exp(-E/RT^)  +  exp  (-E/RT.  )  ] ,  v;hich  in  fact 
averages  the  rate  constant  as  opposed  to  the  temperature  might  be  used. 
(E  =  activation  energy  of  the  rate  controlling  step,)   This  iteration 
across  a  time  step  would  not  significantly  alter  the  results,  and  so 
this  was  not  done  in  this  study. 

Tlie  term  /  h.xj.  in  the  energy  equation  was  treated  similarly. 
Since  the  temperature  was  assigned  its  value  at  the  beginning  of  a 
main  time  step  and  regarded  as  constant  over  the  main  tim.e  step,  the 
h.  (i=l,2, N)  are  also  constant  over  this  main  step  and  so. 
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So  the  result  of  the  previous  integration  was  used  here  instead  of 
resorting  to  a  separate  integration  of  the  term  )_  h.v;.  .   Again,  a 
mean  value  for  the  h.  can  be  used  for  improved  accuracy. 

4.3   Chemical  Kinetics 


The  conservation  equations  contain  terms  involving  the  rate  of 
production  of  species  due  to  chemical  reactions.   These  terms  must  be 
evaluated  before  the  equations  can  be  solved.   The  evaluation  of  these 
terms  requires  a  knowledge  of  the  reaction  mechanism  and  the  rate 
constant  parameters  for  the  different  elementary  reaction  steps. 

Rate  constants  for  the  elementary  reactions  occurring  in 
the  C-O-H-N  system  have  been  tabulated  by  several  authors.   However, 
there  is  considerable  disagreement  over  the  values  of  the  rate 
constants  from  author  to  author.   Again,  there  is  disagreement  over 
the  reaction  mechanism  proposed  for  the  oxidation  of  methane.   In 
this  analysis  a  total  of  16  species"  were  considered.   It  was  found 
that  it  is  possible  to  write  a  total  of  84  (42  forward  and  42  reverse) 
bimolecular  elementary  reaction  steps  for  these  16  species.   In  addition 
to  these,  several  third  order  reactions  involving  a  third  body  are  also 
possible.   The  kinetic  data  for  all  these  reactions  are  not  available; 
moreover  it  is  not  necessary  to  consider  all  the  reactions  for  which 
data  are  available.   A  detailed  analysis  of  all  the  bim.olecular 
reactions  for  which  data  are  available  was  made  to  determine  the 
contribution  of  each  reaction  towards  the  total  rate  of  production 
of  each  of  the  species.   This  was  done  over  the  temperature  range 
600-2400 °K.   On  the  basis  of  this,  reactions  which  contributed  less 


*C0    CO,  N    HO,  H,  OH,  0,  KO ,  NO^ ,  N^O,  N,  0^ ,  H^ ,  CH^, 
CHO,  CH^, 
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than  0.1%  to  the  total  rate  of  production  of  any  of  the  species  over 
the  entire  temperature  range  were  eliminated  from  further  consideration. 
Tl-i.e  contribution  of  each  reaction  from  an  energy  standpoint  was  not 
tested;  it  being  assumed  that  the  reactions  which  contribute  less 
than  0.1%  to  the  total  rate  of  production  of  any  of  the  species 
would  not  cause  any  significant  contribution  to  the  energy  equation 
also,  unless  their  heats  of  reaction  are  at  least  an  order  of 
magnitude  more  than  those  of  the  reactions  which  were  not  eliminated 
from  consideration.   A  similar  analysis  was  also  done  for  the  third 
order  reactions.   Finally,  a  mechanism  consisting  of  25  elementary 
reaction  steps  (Table  4.1)  was  used  for  the  analysis.   Tne  first 
14  reactions  are  the  ones  which  Seery  and  Bowman  [14]  have  suggested 
in  their  proposed  scheme  for  the  oxidation  of  methane.   The  remaining 
reactions  are  concerned  with  the  N-0  system. 

The  kinetic  data  used  in  this  analysis  have  by  no  means  been 
firmly  established,  and  suggested  error  limits  foi"  some  of  the  rate 
constants  are  as  high  as  25%.   The  error  is  greater  the  higher  the 
temperature,  since  most  of  the  experimental  data  are  gathered  in 
the  temperature  range  300-1000  degrees  K  and  involve  extrapolation 
to  higher  temperatures.   Also  the  oxidation  mechanism  of  methane 
has  not  been  firmly  establislied  and  m.uch  work  remains  to  be  done  on 
the  effects  of  pressure,  tempei-ature ,  mixture  composition,  and  the 
presence  of  inerts  on  the  oxidation  mechanism.   In  the  absence  of 
better  data  and  a  more  definite  reaction  mechanism,  the  data  referred 
to  above  were  used.  The   method  developed  however  is  general,  and  as 
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Table   4.1 


Reactions    considered   in   this    analysis 


k      =  AT     exp(-E/RT) 


Forward   rate    constant 


Reaction 

CH,    +  M  = 
4 

=   CH     +   H  +  M 

OH  +  CHO 

=    CO  +   H2O 

°2   +  ^^3 

=    HO  +  CHO 

CHO  +  M  = 

=   CO  +  H  +  M 

0  +   CH,    -- 
4 

=    OH  +   CH 

CH,    +   H  = 
4 

=   H2   +  CH^ 

OH  +  CH, 
4 

=   H^O  +  CH^ 

OH  +   CO   = 

=    CO2   +  H 

H  +   0^    = 

OH  +  0 

0  +  H^    = 

H  +  OH 

0  +  H^O   = 

=    OH  +  OH 

H  +  H^O   = 

=    H^   +   OH 

M  +   H^O  = 

=    H  +   OH  +  M 

0  +   CH      = 

=    HCO   4-   H^ 

N^   +  0   = 

NO  +  N 

NO  +   0  = 

O2   +  N 

N^   +   0^   = 

=    NO   +   NO 

A 


E/R 


Ref. 


NO  +  0     =   0  +  NO 

M  +  NO      =    0  +  NO  +  M 

0   +   NO   =    NO  +   NO 


NO  +  NO^ 


N^O  +   0^ 


NO  +  N      +   O2    =    N^O  +   NO2 
OH  +  M  =    0  +  H  +  M 
H  -f   NO   =   OH  +  N 
N     +  OH  =    H  +  NO 


1.5x10 
1.0x10 
2.0x10 
2.0x10 
1.7x10 
6 . 3x10 
2.8x10 
5.6x10 


19 
14 
10 
13 
13 
13 
13 


11 


2.24x10 


14 


1.74x10 

14 
8.4x10 


13 


1.0x10 
2.0x10 
1.0x10 


14 
22 
14 


1.36x10 
1.55x10' 


14 


4.55x10 

12 
1.0x10 


24 


1 . 1x10 


16 


3.24x10 

14 
6.0x10^^ 


12 


2.30x10 

21 
1.0x10 


16 


5.0x10 
4.0x10 


11 
14 


0 

50300 

0 

0 

0 

0 

0.50 

14400 

0 

4380 

0 

6350 

0 

2500 

0 

543 

0 

8450 

0 

4750 

0 

9120 

0 

10200 

•1.0 

61500 

0 

0 

0 

37700 

1.00 

19450 

■2.50 

64300 

0 

22700 

0 

32500 

0 

32140 

■1.50 

4950 

0 

17885 

•1.50 

50500 

0.50 

2500 

0 

8000 

14 
14 
14 
14 
14 
14 
14 
15 
15 
15 
14 
14 
16 
17 

15 

15 
15 
15 
15 
16 
16 
16 
16 
16 
16 


Units:   E/R  in  degrees  K 
rate  constants  in 
where  N  =  order  of  the  reaction 


(cc/g.mole;    /sec. 
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Table  4.1  (Continued) 


Reactions  which  were  not  considered  because  of  negligible  contribution 

Forward  rate  constant 


Reaction 


n 


E/R 


Ref . 


N  +  NO  =  NO  +  NO 

N^  +  0^  =  0  +  N^O 

NO  +  H  =  HO  +  NO 

OH  +  NO  =  H  +  NO 

0  +  H  =  OH  +  OH 

0  +N  0=  NO  +  N 

N  +  M  =  N  +  N  +  M 

NO   +  NO   =  NO  +  NO  +  0 

NO  +  N  +  0   =  NO  +  NO 

0  +  0  +  M  =  0  +  M 

CO  +  NO^  =  NO  +  CO 

NO  +  NO  =  NO  +  N 

NO  +  M  =  0  +  N  +  M 


4.7x10 

4.0x10 

2.0x10' 

1.0x10 

1.0x10 

3.7x10 


14 
19 


14 
14 
14 


2.97x10 

13 
1.0x10 


21 


2 . 3x10 


16 


1.00x10 
S-OxlO-*-^ 


14 


NO. 


M  =  N  +  0, 


M 


1.00x10 
2.2  7x10 
6.00x10 


10 
17 

14 


0 

41650 

0 

46600 

0 

9000 

0 

600 

0 

35000 

■0.5 

20900 

•1.5 

112500 

0.5 

15500 

0 

17885 

0 

0 

0 

13900 

0 

49300 

■0.50 

74900 

■1.50 

52600 

16 
16 
16 
16 
16 
16 
16 
16 
16 
18 
19 
15 
15 
19 


N^  +  0^ 


NO, 


2.70x10 


14 


-1.00 


60600 


19 


Reactions  for  which  kinetic  data  were  not  available 


CO  +  N^  =  CO  +  NO 
CO   +  OH  =  0   f  CHO 
CO  n-   N  =  CO  +  NO 
CO  +  OH  =  0  +  CHO 


CO  +  CH,  =  CH.,  •+-  CHO 
N  +  NO  =  NO  +   N 
HO  +  N  =  NO  •+  H^ 


CO   +  H  =  0  +  CHO 
CO   +  0  =  CO  +  0 
CO^  +  H^  =  CO  +  H2O 
CO  •+-  H  =  H  +  CHO 
N^  +  H^O  =  N^O  +  H^ 
H^O  +  0  =  0^  -+-  H^ 


Most  of  these  reactions  are  not  of  much  significance  from  a 
practical  standpoint. 
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more  is  known  about  the  oxidation  mechanism  and  the  rate  constants, 
the  new  information  can  be  readily  incorporated  to  obtain  more 
realistic  solutions.   In  the  computer  program  developed  any  number 
of  species  and  elementary  reaction  steps  may  be  considered,  the  only 
limi.tation  being  the  com.puter  time  and  storage  available. 

4.4  Uncertainty  in  the  Definition  of  Temperature 

It  is  now  generally  realized  that  there  does  exist  a  distinct 
possibility  that  the  various  chemical  species  in  a  flame  zone  may  not 
be  in  equilibrium  with  respect  to  their  various  degrees  of  freedom. 
Thus  the  concept  of  a  'temperature'  to  define  the  thermodynamic 
state  of  the  gases  in  the  flame  zone  is  rendered  difficult.   This 
problem  is  briefly  discussed  here;  details  may  be  found  in  the 
appropriate  literature. 

llie  statistical  definition  of  temperature  is  based  on  the 

Maxwell-Boltzmann  equilibrium  energy  distribution  function.   It  is 

defined  as  a  parameter  characteristic  of  a  system  at  equilibrium 

and  it  specifies  the  distribution  of  energy  between  the  particles 

comprising  the  system.   Normally  the  definition  is  based  on  the 

distribution  of  the  translational  kinetic  energy.   Since  polyatomic 

molecules  have  in  general  other  modes  of  energy  storage,  it  is 

possible  to  analogously  define  various  temperatures,  e.g.,  T 

T    ,  T  ^   .   For  a  system  at  comnlete  equilibrium  all  these 
rot'   elec 

definitions  give  a  unique  value  of  temperature,  the  various  degrees 
of  freedom  being  at  equilibrium  with  each  other. 
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The  kinetic  theory  definition  of  temperature  for  a  system 
composed  of  different  chemical  species  is  expressed  in  terms  of  the 
mean  particle  energy,  averaged  over  all  the  species  [3,  p.  455]: 


f  «  4  J^  n.  (i  ».V^^) 


Temperature  for  the  other  degrees  of  freedom  m.ay  be  defined 
analogously.   For  a  multi-component  system  not  at  equilibrium,  it  is 
possible  that  not  only  are  the  various  degrees  of  freedom  of  a 
particular  species  not  in  equilibrium  with  each  other,  but  that  there 
is  non-equilibrium  between  the  different  chemical  species,  i.e.,  it 
may  not  be  possible  to  define  a  unique  translational  (or  rotational 
or  vibrational)  temperature  for  the  entire  system. 

In  a  laminar  premixed  flame,  the  temperature  of  the  gas  is 
raised  extremely  rapidly  (more  so  in  detonations  and  shock  waves). 
Ivhen  this  happens  the  different  degrees  of  freedom  would  in  general 
come  to  equilibrium  at  different  rates.   Translational  equilibrium 
is  attained  the  fastest — in  about  five  mean  free  paths,  as  demonstrated 
by  shock-wave  experiments  in  gases  with  no  internal  degrees  of  freedom 
(Ne,  Ar,  Kr) .   The  temperature  of  primary  concern  here  is  the  trans- 
lational one,  since  this  is  used  in  the  equation  of  state  for 
transport  property  calculations  and  in  the  evaluation  of  reaction 
rate  constants.   Since  the  assumption  of  a  continuum  has  already 
been  made,  i.e.,  changes  in  macroscopic  properties  occur  over 
distances,  large  com.pared  to  a  mean  free  path,  the  concept  of  a 
local  translational  temoerature  should  be  valid. 
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The  rotational  degrees  of  freedom  in  diatomic  and  polyatomic 
molecules  also  come  to  equilibrium  at  a  rapid  rate  though  not  as  fast 
as  the  translational  mode.   It  is  the  vibrational  degree (s)  of  freedom 
which  poses  problems.   Equilibration  or  relaxation  times  for  the 
vibrational  degrees  of  freedom  are  in  general  much  larger  than  for 
the  translational  or  rotational  degrees.   The  reason  for  this  is  that 
energy  differences  between  quantized  vibrational  energy  levels  are 
much  greater  than  for  the  translational  or  rotational  levels.   As 
a  result,  when  the  temperature  is  being  raised  rapidly,  the  numier 
of  collisions  needed  to  excite  the  next  higher  vibrational  energy 
level  may  be  significant,  of  the  order  of  5000,  compared  to  50  for 
translational  and  rotational  equilibrium.   For  triatoraic  molecules 
(NO,  CO  )  it  may  be  more,  about  2500-50,000.   Thus  the  translational 
and  rotational  temperatures  of  a  gas  being  heated  rapidly  follox-/ 
their  equilibrium  values  closely  V7hile  the  vibrational  temperature 
may  still  be  its  cold-gas  value.   In  such  cases  a  strict  unique 
definition  of  temperature  is  rendered  impossible.   It  may  be  neces- 
sary to  define  two  different  tem.peratures — one  translational  and  the 
other  vibrational. 

Thus  the  definition  of  temperature  is  linked  to  the  requirement 
of  equilibrium.   If  there  is  non-equilibrium  between  the  different 
chemical  species  or  the  different  degrees  of  freedom  of  a  particular 
species,  the  Maxwell-Boltzm.ann  equilibrium  distribution  function 
ceases  to  hold,  and  so,  strictly  speaking,  a  temperature  cannot  be 
defined  for  the  system.   To  use  the  concept  of  temperature  for  a 
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non-equilibrium  system,  the  assumption  generally  made  is  that  even 
though  the  system  as  a  whole  is  not  at  equilibrium,  local  equilibrium 
exists  at  any  particular  point  in  the  system.   This  boils  do^^n  to 
assuming  that  in  any  elementary  volume  (whose  dimensions  are  'small' 
compared  to  those  of  the  system)  the  particles  of  the  different  species 
are  in  equilibrium  with  respect  to  each  other  and  their  individual 
degrees  of  freedom,  so  that  the  Maxwell-Boltzmann  energy  distribution 
may  be  used  to  define  a  local  temperature.   This  is  how  the  use  of 
temperature  in  flame  analysis  is  justified. 

In  flames  where  the  transfer  of  energy  between  the  vibrational 
and  trans lational  energy  modes  is  extremely  slow,  the  need  for  separate 
definitions  of  translational  and  vibrational  temperatures  may  destroy 
the  concept  of  a  local  temperature  also.   Such  cases  have  been  shovm 
to  exist  by  spectroscopic  m.eans  [20].   Ttie  chemical  reaction  rates 
and  the  transport  properties  are  evaluated  at  values  of  temiperature 
based  on  a  static  system  where  local  equilibrium  has  been  achieved. 
For  the  case  just  mentioned  these  evaluations  may  be  somewhat  different 
due  to  non-equilibrium,  effects;  it  is  not  known  whether  the  difference 
is  significant  in  all  cases.   Som.a  work  on  the  high  temperature 
oxidation  of  hydrogen  and  the  dissociation  of  oxygen  has  been 
reported.   The  results  show  that  the  slowness  of  vibrational  relaxa- 
tion does  change  the  processes  to  a  certain  extent;  the  chemical 
reactions  occur  at  'temperatures'  different  from  their  local 
equilibrium  values. 

For  chemical  reactions  \7hich  are  extremely  rapid  (half-lives 
of  10  microseconds  or  less)  the  classical  assum.ption  of  equipartition 
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of  energy  is  invalid.   Observations  for  some  rapid  reactions  have 
revealed  that  some  of  the  reaction  products  are  formed  almost 
entirely  in  higher  vibrational  levels  (e.g.,  14th  or  15th)  and  the 
transition  to  the  equilibrium  ground  state  occurs  rather  slowly  [21] 
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CHAPTER  V 
l-IATHEtlATICAL  SBIULATION  OF  THE  OTTO  CYCLE 

5.1   Introduction 

The  Otto  cycle  may  be  envisioned  as  a  sequence  of  thermodynamic 
processes,  during  which  the  thermodynamic  properties  and  the  cheirdcal 
composition  of  the  working  fluid  undergo  changes.   The  cycle  starts 
with  the  intake  of  a  fuel-air  mixture  at  state  1  (Fig.  5.1).   The 
mixture  is  compressed  until  state  2,  where  it  is  ignited  with  an 
electric  spark.   Process  2-3  represents  the  combustion  and  3-4  is  the 
expansion  process  at  the  end  of  which  the  gases  are  expelled  from  the 
engine  cylinder  and  a  fresh  charge  drax\m  in. 

The  early  theoretical  analyses  made  several  assum.ptions  to 
simplify  the  problem  and  make  it  solvable  without  a  great  deal  of 
computational  work.   For  example,  it  is  common  to  assume  that 
processes  1-2  and  3-4  are  isentropic  and  that  the  combustion  process 
2-3  and  exhaust  process  4-1  occur  at  constant  volume.   Further 
simplifications  are  obtained  through  the  assumption  that  air  is  the 
working  fluid  throughout  the  cycle  and  its  properties  (specific  heats) 
are  those  of  atmospheric  air  (known  as  the  air  standard  cycle).   Recent 
analyses  have  removed  several  of  these  simplifying  assumptions  and 
modeled  the  processes  occurring  in  the  Otto  cycle  engine  on  a  more 
realistic  basis.   'This  has  largely  been  possible  due  to  the  availability 
of  high-speed  electronic  computing  machines.   A  review  of  the  early 
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Fig.    5.1.      p-V   diagram   for  the   Otto    cycle 
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work  is  given  in  Ref.  1.  Among  some  of  the  recent  work  considering 
the  complex  sequence  of  processes  to  some  detail  are  those  of  Edson 
[2],  Spadacclni  [3],  and  others. 

Edson 's  analysis  was  made  on  the  assumption  of  chemical 
equilihrium  at  every  step  in  the  cycle;  the  main  aim,  however,  was 
to  study  the  influence  of  compression  ratio  and  dissociation  on  the 
thermal  efficiency.   Tnis  assumption  of  chemical  equilibrium  has 
been  made  in  several  analyses  due  to  the  lack  of  chemical  kinetics 
data. 

Spadaccini's  analysis  [3]  has  been  one  of  the  most  detailed 
analyses  made  so  far.   In  this  analysis  the  compression  process  is 
assumed  to  be  isentropic  and  the  specific  heat  is  assumed  to  be  a 
linear  function  of  temperature  for  this  process.   Combustion  is 
assumed  to  occur  at  constant  volume  and  the  composition  and 
thermodynamic  properties  of  the  mixture  after  combustion  are  found 
on  the  basis  of  a  chemical  equilibrium  analysis.   For  the  expansion 
process,  a  non-equilibrium  finite  rate  kinetics  analysis  with 
allowances  for  heat  transfer  is  m.ade  considering  a  total  of  13 
chem.ical  species  to  be  present.   It  will  be  seen  later   that  an 
equilibrium  analysis  of  the  expansion  process  does  not  convey  the 
true  state  of  affairs  with  respect  to  the  concentrations  of  the 
species  at  exhaust.   Since  the  present  interest  is  to  study  the 
exliaust  emissions,  a  non-equilibrium  analysis  is  a  must  to  obtain  a 
realistic  estlm.ate  of  the  emission  levels.   Other  work  in  this  area 
has  been  referred  to  at  the  appropriate  places  in  this  chapter. 
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5.2   Present  Work 

Here  again  the  compression  process  has  been  assumed  to  be 
isentropic.   The  combustion  process  is  initiated  sometime  before  top 
dead  center  is  reached  and  completed  shortly  after  top  dead  center. 
Instead  of  assuming  the  combustion  process  to  be  instantaneous  as 
done  in  many  of  the  previous  analyses,  here  the  process  has  been 
assumed  to  proceed  at  a  finite  rate.   Using  an  assumed  burning  rate 
law  which  gives  the  mass  fraction  of  the  total  fuel-air  mixture  burnt 
as  a  function  of  the  crank  angle,  a  step-by-step  analysis  of  the 
combustion  process  has  been  made  to  obtain  the  state  properties  and 
the  chemical  composition  of  the  working  fluid  at  the  point  where 
combustion  is  complete.   After  this,  for  the  expansion  process  analysis, 
the  same  reaction  mechanism  as  used  by  Spadaccini  has  been  used  here; 
however,  different  mathematical  techniques  have  been  used  to  solve  the 
governing  equations. 

The  concentrations  of  the  various  species  at  exhaust  are  thus 
obtained  and  so  this  model  can  be  used  to  study  the  influence  of  several 
parameters  such  as  air-fuel  ratio,  exhaust  gas  recirculation,  compression 
ratio,  speed,  ignition  timing,  etc.,  on  the  exhaust  endssions.   Due  to 
the  approximations  involved  in  the  theoretical  analysis,  and  the  fact 
that  the  processes  occurring  in  the  exhaust  pipe  have  not  been  considered 
in  this  analysis,  agreement  with  the  experimentally  measured  values  may 
not  be  exact.   However,  the  trend  In  the  results  will  definitely  show 
up  in  this  theoretical  analysis  and  so  experimental  runs  need  be  made 
to  obtain  results  only  for  a  few  combinations  of  the  different 
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parameters  which  seem  promising  from  the  emissions  viev7point ,  after 
a  preliminary  detailed  theoretical  analysis. 

The  only  empirical  input  here  is  the  burning  rate  law.   The 
validity  of  the  burning  rate  law  can  be  checked  by  comparing  the 
calculated  thermodynamic  properties  (i.e.,  the  pressure  and  the 
temperature)  with  the  experimentally  measured  values.   The  law  may 
then  be  modified,  if  necessary,  to  obtain  better  agreement  with 
experimental  results.   The  assumptions  made  here  agree  qualitatively 
with  the  experimental  results  of  Refs.  4  and  5  and  other  published 
literature;  quantitative  agreement  of  the  mass  fraction  burnt  as  a 
function  of  the  crank  angle  could  not  be  checked  due  to  differences 
in  the  operating  conditions. 

A  rigorous  theoretical  analysis  involving  the  solution  of  the 
general  three-dimensional  time-dependent  conservation  equations  for 
a  flame  propagating  into  a  premixed  fuel-air  mixture  would  yield  the 
mass  fraction  burnt  as  a  function  of  the  crank  angle.   However,  with 
present  day  computational  techniques,  this  type  of  an  analysis  of 
the  problemi  is  almost  impossible.   Moreover,  the  analysis  is  beset 
by  difficulties  due  to  the  limited  knowledge  of  the  com.plex  nature 
of  the  processes  occurring  in  the  engine:   the  reaction  mechanism 
aad  the  rate  constant  values,  degree  of  turbulence  in  the  combustion 
chamber,  problems  arising  from  the  fuel-air  mixture  not  being  homo- 
geneous and  the  fuel  not  being  completely  vaporized,  heat  transfer 
at  the  combustion  chamber  walls  and  the  associated  quenching  effects 
on  the  flame,  together  with  a  host  of  other  difficulties. 
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The  effect  of  the  injection  of  alcohol-water  mixtures  on  the 
emissions  has  been  studied  using  the  model.   In  general,  many  of  the 
control  techniques  which  have  been  suggested  involve  som.e  treatment 
of  the  exhaust  gases.   The  injection  of  water-alcohol  mixtures  strives 
to  reduce  the  formation  of  the  pollutants  in  the  engine  itself  rather 
than  offer  some  treatment  of  the  exhaust  gases. 

The  hydrocarbon  emissions  result  mainly  from  the  heterogeneous 
processes  occurring  in  the  quench  zone  and  are  influenced  by  several 
other  factors.   Due  to  the  lack  of  knowledge  of  the  complex  processes 
occurring,  it  was  not  possible  to  predict  the  hydrocarbon  emissions 
using  this  model.   However,  it  should  be  remembered  that  control 
techniques  v/hich  result  in  the  reduction  of  carbon  monoxide  also 
generally  result  in  the  reduction  of  hydrocarbons.   So  the  latter  need 
not  be  studied  separately,  but  can  be  determined  by  subsequent  experiments 
after  the  preliminary  theoretical  analysis. 

5.3  Major  Assumptions  for  Analysis 

Before  describing  the  mathematical  model  of  the  Otto  cycle, 
a  discussion  of  the  assumptions  made  in  the  analysis  is  in  order. 
In  this  section  the  assumptions  relating  to  the  entire  cycle  are 
mentioned.   Assumptions  for  the  specific  processes  are  m.entioned  in 
the  appropriate  places  in  this  presentation. 

(a)   The  fuel  is  assumed  to  be  completely  vaporized  and  mixed  homo- 
geiieously  with  air.   If  the  fuel-air  mixture  is  not  homogeneous,  it 
will  be  locally  fuel-rich  at  some  points  and  locally  fuel-lean  at 
others,  though  it  m.ay  be  overall  stoichiometric.   Compared  to 
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stoichiometric  combustion,  the  fuel-rich  parts  would  yield  larger 
amounts  of  unburnt  hydrocarbons  and  carbon  monoxide  while  the 
fuel-lean  parts  would  yield  a  larger  amount  of  the  oxides  of 
nitrogen.   Also  the  energy  release  would  be  less  since  a  lesser 
amount  of  the  fuel  is  actually  burnt.   The  result  would  be  (as 
compared  to  the  homogeneously  mixed  case) : 

(i)   relatively  low  mean  exhaust  temperature,  and 
(ii)   higher  concentrations  of  carbon  monoxide,  unburnt 
hydrocarbons,  and  oxides  of  nitrogen. 

Similar  conclusions  would  hold  when  the  mixture  was  overall 
fuel- lean  or  fuel-rich,  and  not  homogeneous.   Thus  the  theoretical 
calculations  assuming  perfect  mixing  represent  the  upper  limit  of 
performance  and  the  lower  limit  of  the  emissions. 

(b)  The  gas  mixture  is  assumed  to  be  thermally  perfect.   Nitrogen 
is  a  major  constituent  (>  70%)  of  the  mixture  and  it  closely 
approximates  thermally  perfect  behavior  for  the  conditions 
encountered  in  a  typical  Otto  cycle  engine.   Other  constituents  are 
present  in  relatively  smaller  concentrations,  so  as  not  to  affect  the 
overall  thermally  perfect  behavior.   The  only  constituent  likely  to 
be  considerably  thermally  imperfect  would  be  a  fuel  like  iso-octane. 
However,  since  its  concentration  in  the  mixture  is  small  (mole  fraction 
<  0.02),  the  results  would  be  affected  only  slightly. 

(c)  Reactions  in  the  unburnt  gas  region  are  neglected,  i.e.,  there 
are  no  precombustion  reactions.   The  maximum  tem.perature  in  the 
unburnt  gas  region  is  generally  less  than  1000  degrees  K  for  normal 
combustion,  so  this  is  a  fairly  good  assumption.   Some  reactions 
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(i.e.,  low  temperature  pyrolysis,  slow  oxidation,  cool  flames)  do 
occur  in  this  region  in  practice,  but  the  rates  of  these  reactions 
are  significant  only  at  temperatures  above  about  1400  degrees  K. 

(d)  The  quench  layer  is  assumed  to  occupy  a  negligible  fraction  of 
the  total  volume.   The  quench  layer  is  a  major  source  of  the  hydro- 
carbons and  products  of  incomplete  combustion  present  in  the  engine 
exhaust.   This  assumption  implies  that  only  the  homogeneous  chemical 
and  physical  processes  occurring  in  the  bulk  of  the  gases  are 
considered  and  the  heterogeneous  processes  occurring  near  the  combustion 
chamber  walls  are  neglected. 

(e)  There  is  no  heat  flow  or  mixing  between  adjacent  elements  of  gas 
during  combustion,  i.e.,  transport  phenomena  of  viscosity,  thermal 
conductivity  (also  radiation  and  convection)  are  neglected.   Transport 
phenomena  are  of  importance  only  in  the  flame  element  where  appreciable 
gradients  of  temperature  and  species  concentrations  exist,  and  these 
have  been  considered  in  the  analysis  of  the  flame  zone  presented 
earlier  in  this  study.   In  the  unbumt  gas  region,  gradients  in  these 
properties  are  generally  absent,  since  this  region  is  fairly  homogeneous. 
In  the  burnt  gas  region,  gradients  in  these  properties  are  present,  but 
are  not  appreciable  enough  to  cause  significant  fluxes  of  heat  and  mass. 
Besides,  in  this  study  a  bulk  analysis  has  been  made  considering  the 
entire  burnt  gas  region  to  be  at  a  uniform  temperature  and  concentration. 

(f)  No  pressure  gradients  exist  in  the  combustion  chamber.   Although 
the  pressure  changes  with  time,  it  is  spatially  uniform  within  the 
combustion  chamber;  any  pressure  gradient  introduced  v7ould  be 
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equalized  at  the  velocity  of  sound  in  the  medium,  which  is  large 
compared  to  the  flame  velocity.   This  would  not  be  true  under 
knocking  conditions. 

(g)   Solid  phases  (e.g.,  solid  carbon)  are  assumed  to  be  absent  in 
the  burnt  products.   This  would  be  a  reasonable  assumption  except 
when  the  mixture  is  extremely  fuel-rich  either  on  an  overall  or  on 
a  local  basis.   These  cases  are  not  considered  here, 
(h)   The  burnt  gas  region  is  assigned  one  single  temperature  and 
composition.   In  practice,  a  temperature  and  concentration  (of  some 
or  all  species)  gradients  are  present;  higher  temperatures  being 
nearer  the  point  of  ignition.   Tlie  gradient  is  neglected  and  the 
uniform  temperature  and  com.position  are  determined  by  an  averaging 
process.   This  is  discussed  further  in  Section  5.6. 

5. A   Theoretical  Analysis 

1.   Calculations  of  conditions  at  beginning  of  compression. 
The  first  step  in  the  cycle  analysis  is  an  estimation  of  the  mass 
fraction  (f),  the  temperature,  pressure  and  composition  of  the 
residual  gases  in  the  cylinder  from  the  previous  cycle.   The  m:ass 
fraction  of  the  residual  gases  is  obtained  readily  using  the  following 
expression  [6]: 


f  =        ^  ■      ..-rr  (1) 


T  +  T 
1    4 


Pi 
r  — 
P 
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4    ^P4^ 


where  r  =  compression  ratio,  and  the  other  symbols  have  their  usual 


meanings , 
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The  temperature  and  composition  of  the  residual  gas  are 
estimated  on  the  basis  of  previous  work  reported  in  this  area  and 
the  pressure  is  assumed  to  be  the  same  as  the  intake  pressure,  namely, 
atmospheric.   The  value  of  f  varies  from  0.01  to  0.1  (generally 
0.05-0.06),  so  errors  in  the  estimation  of  the  temperature  and  chemical 
composition  will  not  seriously  affect  the  remainder  of  the  analysis. 

The  conditions  at  the  beginning  of  the  compression  stroke  are 
obtained  by  assuming  that  the  residual  gases  and  the  incoming  fuel-air 
mixture  mix  at  a  constant  pressure  without  heat  transfer  to  the 
surroundings  (i.e.,  at  constant  total  enthalpy)  and  without  chemical 
reactions.   The  result  is  a  homogeneous  m.ixture  at  one  uniform 
pressure  and  temperature  whose  enthalpy  is  equal  to  the  sum  of  the 
enthalpies  of  the  residual  gases  and  the  fuel-air  mixture.   The 
composition  of  the  mixture  is  given  by: 
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where  the  subscripts  'ex'  and  'mix'  refer  to  the  residual  gases  and 
the  incoming  fuel-air  mixture  respectively. 

A  Newton-Raphson  iteration  scheme  is  used  to  calculate  the 
resultant  tem.perature  T  based  on  the  condition  Ah  =  0.   The  relevant 
equations  are: 


h  =  y  x.h. 


(3) 


h(T^,p^)    =  h(Tj^,pp   +  AT|^ 


(4) 


(T     p    )    are    the  kth   estimates   of   temperature    and  pressure   respectively, 
Ti\e  partial   derivative    is   evaluated   at    (T     p,  )    and   is    given  by: 
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The  iteration  is  carried  out  until  the  temperature  converges  to  within 
a  predetermined  accuracy.   This  is  checked  by  the  condition: 


|AT  /T  !  ^  0.0001 


(8) 


The  first  estimate  for  T  is  obtained  from: 
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Tlie  thermodynamic  properties  used  in  the  above  equations  and  later  on 
in  the  analysis  are  all  obtained  from  the  JANAF  Tables. 

2.   The  compression  process.   Having  found  the  thermodynamic 
properties  at  the  beginning  of  the  compression  process  (BDC)  the 
conditions  at  ignition  are  determined  by  assuming  isentropic 
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compression  without  chemical  reactions  from  BDC  to  the  point  of 
ignition.   The  following  iteration  scheme  is  used: 

s=yx.s.  -R  y  x.lnx.    -  R  y  x.^nP    (P  is  in  atm. )     (10) 
'-ii    u^xi    u^i 
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The  partial  derivatives  are  evaluated  at  (T  p  )  and  are; 
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Tlie  first  estimates  of  (T  ,p^)  are  based  on  the  relation  pv  =  constant, 
with  a  suitable  mean  value  for  Y;  a  value  of  1.35  was  used  here. 
^■^(T^jP-,)  is  the  volume  corresponding  to  the  crank  angle  at  ignition 
and  is  obtained  from  the  following  equation  expressing  the  combustion 
chamber  volume  as  a  function  of  the  crank  angle  6: 
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(G  is  measured  from  TDC) 
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where   £  =  crank  radius 

L  =  connecting  rod  length 
PD  =  piston  displacement 
CR  =  compression  ratio 

Equations  (11)  and  (12)  are  solved  simultaneously  for  AT^^,  Ap^  and 
the  next  estimates  are  obtained: 


T,  ^1  =  T,  +  AT,  (18) 

k+1    k     k 


Pk+1  =  Pk  +  ^Pk  ^^^> 


The  process  is  repeated  until  the  temperature  and  the  pressure  converge 
to  within  a  predetermined  accuracy. 

The  above  calculations  give  values  of  the  therm.odynamic 
properties  v;hich  are  sufficiently  accurate  for  the  purpose  here. 
A  more  detailed  analysis  can  be  made  on  the  basis  of  the  wave  motion 
in  the  gases  during  the  compression  process.   As  the  piston  begins  to 
move  in,  it  sets  up  a  series  of  compression  waves  which  raise  the 
temperature  and  the  pressure  of  the  gas.   A  convenient  way  of  analyzing 
the  wave  motion  and  determining  the  resulting  change  in  the  thermo- 
dynamic properties  is  the  well-known  method  of  characteristics.   Details 
of  the  method  may  be  found  in  Ref.  7.   However,  the  amount  of  com.puta- 
tional  work  involved  would  not  warrant  its  use  under  normal  circumstances. 

3.   Combustion  (Process  2-3).   To  obtain  a  realistic  model  of 
the  combustion  process,  it  is  necessary  to  consider  the  chemical 
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reactions  which  produce  the  final  combustion  products.   The  oxidation 
of  hydrocarbons  involves  several  intermediate  species  and  free  radicals. 
High  molecular  weight  hydrocarbons  (such  as  iso-octane)  at  first  undergo 
thermal  decomposition  (pyrolysis)  and  the  products  thus  produced  undergo 
oxidation  later  in  the  process, 

A  detailed  kinetic  mechanism  for  the  oxidation  of  propane  in- 
volving 69  elementary  reaction  steps  between  31  chemical  species  has 
been  suggested  by  Chinitz  and  Baurer  [8].   Mechanisms  for  the  oxidation 
of  methane  have  also  been  suggested.   As  the  hydrocarbon  molecule 
increases  in  complexity,  the  number  of  reaction  steps  and  chemical 
species  involved  increases  rapidly  so  as  to  render  detailed  calculations 
prohibitively  expensive  in  terms  of  computer  running  times. 

A  study  of  the  reaction  times  for  hydrocarbon-air  mixtures 
reveals  that  these  times  are  not  significant  compared  to  the  time 
scale  of  the  processes  occurring  in  an  internal  combustion  engine. 
For  example,  for  an  engine  operating  at  2000  revolutions  per  minute, 
the  time  required  for  one  degree  of  crank  revolution  is  83.3  micro- 
seconds. The   time  required  for  the  completion  of  combustion  is  of 
the  order  of  1000  microseconds  at  a  pressure  of  one  atm.osphere 
(Fig.  5.2)  and  this  time  varies  inversely  with  pressure   so  that  for 


More  specifically,  experiments  indicate  a  correlation  of  the 
type  [9] 

n  =  Vp  '     (for  iso-octane)  (20) 

where   n  =  moles  of  air  in  the  fuel  mixture  consumed  per  unit  time 
V  =  reactor  volume 
p  ~   pressure 

The  pressure  exponent  may  vary  slightly  for  different  fuels. 
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Fig.  5.2.   Chemical  reaction  times  for  hydrogen-air  and  hydrocarbon-air 
mixtures  at  one  atmosphere  pressure  and  unity  equivalence 
ratio  (Ref.  8) 
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the  conditions  existing  in  a  typical  internal  conibustion  engine 
(pressure  =  10-80  atm. ) ,  the  reaction  time  and  the  time  required  for 
one  degree  of  crank  revolution  would  be  of  the  same  order  of  magnitude. 

Thus  it  is  reasonable  to  assume  that  the  combustion  process 
produces  species  which  are  in  chemical  equilibrium.   It  will  be  seen 
later  that  this  assumption  is  violated  by  nitric  oxide  (NO),  and  so 
the  formation  of  nitric  oxide  is  considered  separately  on  a  rate- 
dependent  basis.   Since  the  concentration  of  nitric  oxide  is  small, 
the  calculation  of  the  thermodynamic  properties  and  concentrations  of 
the  remainder  of  the  species  on  the  basis  of  chemical  equilibrium  will 
be  affected  only  slightly  by  the  fact  that  nitric  oxide  is  not  at 
its  equilibrium  concentration. 

Edson  [10]  has  developed  a  mathematical  model  for  combustion 
which  is  suitable  for  considering  the  processes  occurring  in  a  spark 
ignition  engine.   The  combustion  process  is  considered  to  be  m.ade  up 
of  three  steps:   a  combustion  step,  a  piston  movement,  and  a  mixing 
operation.   The  first  step  starts  at  a  given  crank  angle  and  the 
three  steps  are  repeated  one  after  another  until  all  the  fuel-air 
mixture  is  consumed. 

a.   Combustion  step:   Tliis  consists  of  the  constant  pressure 
adiabatic  (hence  isenthalpic)  burning  of  a  certain  mass  of  the 
fuel-air  mixture,  referred  to  as  the  flame  element  (Fig.  5.3).   The 
burning  is  followed  by  the  compression  of  the  system  to  the  total 
energy  and  volume  at  the  beginning  of  the  combustion  step.   The  mass 
fraction  of  the  mixture  burnt  is  obtained  by  means  of  an  empirical 


94 


95 


burning  rate  law.   For  the  case  under  study  here,  with  ignition  at 
30  degrees  before  top  dead  center,  combustion  is  assumed  to  be  spread 
over  40  degrees  of  crank  revolution  and  the  mass  fraction  burnt  (f) 
as  a  function  of  the  crank  angle  6  (measured  from  the  time  of  ignition) 
is  (Fig.  5.4): 

f  =  0.00249^  ,  0  <  e  ^  10 

f  =  0.0488  -  0.24  ,  11  ^  e  1  23  (21) 

f  =  1  -  484  exp (-0.3539)  ,   24  ^  9  <  40 

In  the  first  stages,  the  mass  fraction  burnt  per  degree  of  crank 
revolution  gradually  increases  due  to  increasing  flame  radius  and 
velocity;  this  is  followed  by  a  stage  wherein  the  mass  fraction  burnt 
per  degree  is  constant.   Finally,  it  decays  exponentially  as  the 
flame  is  gradually  quenched  out  towards  the  end  of  the  combustion 
process.   For  coraputational  purposes,  the  combustion  process  is 
considered  to  be  complete  when  99.99%  of  the  fuel-air  mixture  is 
consumed. 

The  chemical  composition  and  the  temperature  of  the  flam.e 
element  are  assumed  to  correspond  to  chemical  equilibrium  conditions 
in  accordance  with  the  discussion  above,  and  are  obtained  by  an 
iterative  process  similar  to  that  given  by  equations  (3)  -  (9).   The 
equilibrium  composition  is  calculated  by  the  method  given  in  Appendix  I 
and  the  condition  of  constant  enthalpy  is  used.   If  it  is  desired  to 
include  the  heat  transfer,  the  change  in  enthalpy  should  be  equated 
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Fig.    5.4.      Mass    fraction   of   fuel~air  mixture   burnt    as    a   func- 
tion  of   crank   angle 
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to  the  heat  transferred  from  the  flame  element.   The  problem  here  is 
obtaining  a  reliable  value  of  the  heat  transfer  which  is  a  complex 
function  of  the  combustion  chamber  geometry  and  the  various  processes 
occurring.   Since  the  combustion  process  is  relatively  rapid  and  the 
surface  area  available  for  heat  transfer  is  small,  the  heat  transfer 
is  generally  neglected.   More  serious  is  the  neglect  of  heat  transfer 
from  the  burnt  and  unbumt  gas  fractions,  since  these  have  much  larger 
surface  areas  exposed  to  the  relatively  cold  combustion  chamber  walls 
than  does  the  flame  element.   Mass  and  heat  transfer  (except  possibly 
by  radiation  and  a  little  by  convection)  between  the  unburnt  and  burnt 
fractions  is  negligible  since  the  two  are  separated  by  the  flame 
element  which  has  vanishing  gradients  of  mass  concentrations  and 
temperature  at  its  extremities. 

The  combustion  step  occurs  at  a  fixed  position  of  the  piston, 
i.e.,  at  constant  volume  and  hence  also  at  constant  internal  energy 
in  the  absence  of  heat  transfer.   Therefore,  after  considering  the 
burning  of  the  flame  element,  the  combustion  step  is  completed  by 
restoring  the  system  by  an  isentropic  process  to  the  total  energy 
and  volume  before  the  combustion  step  v;as  started.   The  pressure  and 
temperature  at  the  end  of  the  combustion  step  were  determined  using 
an  iterative  process: 


s(T^,p')  =  s(T.^^,p^)  +  AT.^  ^1     +  Ap^  -^  (22) 


u(T',p')  =  u(T.^,p^)  -l-AT.^ 


f3sl 

[btJ 

P,2 

fas] 
l^pj 

fSu^ 

+  Ap^ 

'3u 
l^pj 

(23) 


P>3 
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The  partial  derivatives  are  as  before  evaluated  at  (T   ,p  )  and  are 


given  by  equations  (13),  (14)  and: 


du 


dT 


=  C  -  -^  (1-a) 
p   w 


(24) 


Su    a  +  3 


l^Pj 


(25) 


where   a 


9logn 


aiosT 


and   e  = 


9logn 


,3logP. 


Also, 


u  =  h  -  R  T/W 
u 


(26) 


In  equations  (22)  and  (23)  the  primed  quantities  refer  to  the  values  at 
the  end  of  the  coirbustion  step  and  the  subscript  j  equals  f,  b,  or  u 
corresponding  to  the  flame  element,  the  burnt  and  the  unbumt  fractions 
respectively.   Writing  equations  (22)  and  (23)  for  each  one  of  these 
fractions  and  expressing  in  matrix  form: 


ds 


dT 


w. 


P,f 
0 

0 


ds 


9T 


P,f 


P,b 
0 

(du 


3 
8T 


P,u 


w 


p,b 


u   BT 


3pJ 


^9s^ 


ds 
3pJ 


T,f 


T.b 


T,u 


P,u  J   -^  ^  ^^ 


rp 


AT 


fk 


AT 


bk 


AT 


uk 


^Pk 


s(T^,p')-s(T^j^,p^) 


s(T^,p')-s(T^j^,Pj^) 


^(■^^'P'^-^^uk'Pk^ 


U  -  y  x..'.u.(T.,  ,p,  ) 


(27) 
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v;here  w.    denotes   the  mass   of   the   appropriate   fraction;   U  is   the 
total   internal   energy  before   the   combustion   step   started,    and   is 
given  by: 


^  =  I  w.u.  (28) 

■      32 
J 


The  set  of  equations  (27)  was  solved  by  Gauss-Jordan  elimination  with 
maximum  pivoting  to  obtain  the  values  of  the  correction  terms  AT 
(j=b,f,u)  and  Ap   in  V7hat  is  essentially  a  Newton-Raphson  iteration 
process  for  a  set  of  simultaneous  algebraic  equations.   The  next 
estimates  were  thus  obtained: 


^J(k+l)  =  V-"''V    (J=^'f'-)  ^29) 


Pk+l  =  Pk  +  ^Pk  (^°> 


and  the  process  was  repeated  until  the  temperatures  and  the  pressure 
all  converged  to  within  a  predetermined  accuracy. 

b.   The  piston  movement  step  isentropically  compresses  the  gases 
to  a  new  volume  V  which  is  determined  by  the  dimensions  of  the  piston 
and  combustion  chamber  and  is  given  by  equation  (17).   The  appropriate 
equations  are  given  below  and  these  are  solved  by  a  method  simdlar  to 
that  used  in  the  corrbustlon  step. 


s(T'.p')  =  s(T.^,P^^)+AT.^  [Il]   _+Ap^ 


P,J 


9s 


op 


(31) 
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V  =  Vm,p')  =  V(T.^,p^)+AT.^  [|I_ 


+  Ap,, 


''sv^ 


P,J 


19PJ 


(j=b,f,u)  (32) 


T,j 


The  partial  derivatives  are  given  by  equations  (13),  (14),  (15),  and  (16), 
Finally, 


di 


3T 


P,f 
0 


(ds 


l3Tj 


p,b 


9s 


3T 


p,u 


f3vl 


w. 


p,f 


9T 


p,b 


ds 


dT 


l^pj 

8pJ 

Ss' 
3p. 

w . 


T,f 


T,b 


T,u 


P,u  J 


'^^fk 

-^^k 

= 

AT  , 
uk 

Ad, 
'k 

s(T^,p')-s(T^^,Pj^) 


s(t;,p')-s(Tj^,^,P^) 


^(T^'P'^-^'^uk'V 


V  -  I   w  v(T    p^) 
3 


(33) 


^j(k+l)  =  ^jk  +  ^"jk    (J=b>f>-) 


(34) 


Pk+1  =  Pk  +  ^Pk  ^^^^ 


(Trie  notation  is  similar  to  that  used  before.) 

c.   Tlie  mixing  operation  consists  of  two  steps.   First,  the 
burnt  gases  and  the  flame  element  are  mixed  together  xjithout  chemical 
reactions  and  the  resulting  burnt  gas  composition  is  obtained  using  a 
relation  similar  to  equation  (2),   To  obtain  the  resulting  pressure 
and  temperature,  two  different  models  were  considered:   constant 
pressure  and  constant  volume.  T\\e.    final  results  were  found  not  to 
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differ  significantly  (less  than  0.02%  difference  in  species  concentra- 
tions and  state  properties)  for  either  mixing  process,  the  differences 
in  the  intermediate  results  not  being  of  interest.   The  second  step 
in  the  mixing  operation  is  the  isentropic  compression  (or  expansion) 
of  the  mixed  burnt  fraction  and  the  unburnt  fraction  to  a  uniform 
system,  pressure  (for  the  constant  volume  mixing  case)  and  the 
specified  total  volume.   The  appropriate  equations  are  given  below  and 
the  same  solution  method  used  as  before. 

First  Step: 

Constant  Volume: 


U  =  w  u^  +  ^^fUf  (36) 

U(T^j^,p)    =    (wj^   +wpu(T^^,p)  (37) 


bk         (9u/3T)       (w^   +  w^)  ^'^^^ 


'^b(k+l)   =   ^bk-^^k  ^39) 


V^bJlVlb 
(V^    +  V^)W 


(40) 


where   W  =  mean  molecular  weight  of  the  mixture  of  burnt  fraction 
and  the  flame  element 

U  =  total  internal  energy  of  the  burnt  fraction  and  the  flame 
element  at  the  beginning  of  the  mixing  process 

R  =  universal  gas  constant 
u 
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Constant  Pressure: 


H  =  w,  h,    +  w  h  (41) 

b    b  ft 


H(Tj^j^,p)   =    (w^   +wph(T^j,,p)  (42) 


H  -  H(T        p) 

AT       =  ~ (43) 

bk        (w^  +  w^)(8h/8T)  ^^^■' 

b  r  p 


^b(k+l)   =   \k-*-^Tbk  <^^) 


R    (w      +  w    )T 

V^        =  V     +     "     \     ^     ^  (45) 

tot  u  pW 


where    H  =  total  enthalpy  of  the  burnt  fraction  and  the  flame 
element  at  the  beginning  of  the  mixing  process 

Second  Step:   Tlie  set  of  equations  (33)  is  used  omitting  the  terms  in- 
volving the  flame  element  and  treating  the  'b'  fraction  to  represent 
the  mixed  'b+f  fraction, 

4.   Expansion  stroke  analysis.   After  combustion  is  complete,  the 
combustion  chamber  contains  a  homogeneous  mixture  of  burnt  products  with 
known  chemical  composition  and  thermodynamic  properties.   As  the  expansion 
process  begins,  the  therm.odynamic  properties  undergo  changes  and  the 
chemical  composition  of  the  system  is  governed  by  finite  rate  kinetics. 
A  set  of  34  elementary  reaction  steps  (17  forward  and  17  reverse)  was 
used  for  the  reaction  kinetics  analysis.   Tliese  are  given  in  Table  5.1 
along  with  the  rate  constant  parameters.   This  reaction  miechanism  is 
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Table  5.1 


Reaction  ■mechanism 


Pveaction 

A 

1, 

0  +  N^   = 

=   NO  +  N 

1.36x10 

2. 

NO  +  0  = 

=   0^  +  N 

l,55xlo' 

3. 

N^  +  0^ 

=   NO  +  NO 

4.55x10' 

4. 

M  +  NO  = 

=  M  +  N  +  0 

2.27x10 

5. 

N^   -f  0^ 

=  N  +  NO 

2.70x10 

6. 

NO  +  0^ 

=   NO^  +  0 

1.00x10 

7. 

NO  +  NO 

=   N  +  NO^ 

1.00x10 

8. 

M  +  NO 

=    0  +  NO  +  M 

1.10x10 

9. 

M  +   NO^ 

=   0     +  N  +  M 

6.00x10 

10. 

CO   +  OH 

=    H  +   CO^ 

5.60x10 

11. 

OH  +   H2 

=    H^O  +   H 

2.19x10 

12. 

OH  +   OH 

=    0  +   HO 

5.75x10 

13. 

0   +  H^    = 

=   H  +  OH 

1.74x10 

14. 

H  +   0^    - 

=    0  +  OH 

2.24x10 

15. 

M  +   0  + 

0=0^+  M 

1.00x10 

16. 

M  +   H  + 

H  =    H     +  M 

7.00x10 

17. 

II  +   H  + 

OH  =   H^O   +  M 

1.50x10 

Forward  rate  constant  parameters 

b       E/R      Ref , 


14 


24 
17 
14 
12 
10 
16 
14 
11 
13 
12 
13 
14 
14 
17 
16 


0 

37700 

1.00 

19450 

2.50 

64300 

0.50 

74900 

•1.00 

60600 

0 

22700 

0 

49  300 

0 

32500 

•1.50 

52600 

0 

543 

0 

2590 

0 

39  3 

0 

4750 

0 

8450 

0 

0 

•1.00 

0 

0 

0 

11 
11 
11 
11 

12 
11 
11 
11 
12 
11 
11 
11 
11 
11 
13 
13 
13 


Forward  rate  constant  =  k   =  AT  exp (-E/RT) 

Reverse  rate  constant  =  (f onward  rate  constant)/ (equilibrium  constant) 

Units:   E/R  in  degrees  K 

rate  constants  m  (cc/g.mole;    /sec. 

V7here  N  =  order  of  the  reaction 
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the  same  as  used  in  Ref.  3.   The  rate  of  production  of  species  'i'  due 
to  chemical  reactions  is  given  by  the  law  of  mass  action,  and  expressed 
in  mathematical  terms,  this  is: 


dt     /,    ij     ij 
J=l 


N    V 


R 


n    n/J   - 


'^  1^1  ' 


N   V 

k    n  n„ 

^J  £=1   ^ 


^j 


(46) 


for  the  M  elementary  reaction  steps  given  by; 


N 


N 


I      V^.R.  =   I   V^.R. 


i=l 


i=l 


(47) 


Equation    (46)   may  be  Xi/ritten   as: 


dn. 


dt 


-  =    fi(^x''^2  ''^3'' 


•  Hj^)        i=l,2,, 


.N 


(48) 


This  is  a  set  of  N  first  order  ordinary  differential  equations.   The 
set  of  equations  was  found  to  be  stiff,  i.e.,  the  time  constants  of 
the  different  equations  differ  by  an  order  of  magnitude  or  more.   As 
such,  standard  integration  procedures  like  Runge-Kutta's  method  are 
not  suitable  since  they  need  prohibitively  long  computer  running 
times  [13].   An  algorithm  developed  by  Gear  [14]  for  the  integration  of 
stiff  equations  was  used  here.   Details  of  the  method  are  given  in 
Appendix  III.   Instead  of  using  a  first  order  starting  procedure  as 
suggested  by  Gear,  a  higher  order  starting  procedure  was  used  here. 
Tills  necessitated  the  calculation  of  higher  order  derivatives  of  n.. 
This  was  done,  using  the  equations  (60a-60f)  given  in  Chapter  IV. 
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Heat  transfer.   The  calculation  of  the  heat  transfer  from 
the  hot  burnt  gases  to  the  cylinder  walls  represents  a  complex 
problem.   In  order  to  make  quick  calculations,  the  problem  would 
have  to  be  simplified  to  such  an  extent  that  the  solution  would  lose 
its  significance.   Overbye  et  al.  [15]  have  made  a  detailed  analysis 
of  the  heat  transfer  in  internal  combustion  engines.   Their  experi- 
mental results  show  that  no  simple  correlation  is  possible. 
Spadaccini  [3]  has  used  an  expression  derived  by  Eichelberg  [16J  for 
heat  transfer  calculations.   In  the  absence  of  a  more  convenient 
relation  this  relation  was  used  here,  although  it  has  been  shoxTO  to 
be  approximate  for  an  Otto  cycle  engine  (the  relation  was  originally 
derived  from  the  data  obtained  for  a  diesel  engine).   The  relation  is: 


4^  =  hA(T   -  T)  (49) 

dt       w 


The  heat  transfer  coefficient  is  given  by: 

h  =  0.0564(V)^/^  (pT)^/^  (50) 

where   V  =  mean  piston  speed  (ft. /sec.) 

A  =  area  exposed  (sq.  ft.) 

T  =  gas  temperature  (degrees  R) 

T  =  wall  temperature  (degrees  R) 
w 

h  =  heat  transfer  coefficient  (BTU/hr.ft.  °R) 
p  =  pressure  (psia) 

The  expansion  stroke  analysis  was  carried  out  in  increments 
of  one  degree  of  crank  revolution  according  to  the  following  scheme: 
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a.  Knowing  the  temperature,  pressure,  and  chemical  composi- 
tion of  the  mixture,  the  reaction  kinetic  equations  were  integrated 
for  a  time  interval  corresponding  to  one  degree  of  crank  revolution, 
assuming  that  the  pressure  and  temperature  remain  constant  over 
this  small  interval. 

b.  This  step  involves  a  piston  movement  corresponding  to 
one  degree  of  crank  revolution.   The  pressure,  temperature,  and 
other  thermodynamic  properties  at  this  volume  [determined  by  the 
dimensions  of  the  engine  according  to  equation  (17)]  were  deter- 
mined using  the  following  iterative  scheme;  this  time  assuming 
that  the  mixture  composition  is  constant  for  this  step: 


u(T    ,p    )  =  u(T,  ,p,  )  +  AT, 
new  -  new       k"^k      k 


9u 


3T 


+  ^P. 


3u 
13PJ 


(51) 


V    =  V(T,  ,p,  )  +  AT   i— ;   +  Ap, 
new     ^  k"^k'     k    \  dT'      ^^■ 


3V 


I'^PJ 


(52) 


where  u  (T    ,p    )  is  obtained  from  the  first  law  of  thermodynairdcs , 
new  ■ new 

knov7ing  the  heat  transfer  and  the  work  done  corresponding  to  this 
piston  movement. 

Tlie  partial  derivatives  are  evaluated  at  (T,  ,p,  )  and  are 
given  by  equations  (15),  (16),  (24)  and  (25).   Also: 


w  =  y  x.w. 
^  11 


T    =  T   +  AT 
k+1    -"k     k 


Pk+1  =  Pk  +  ^h 


(53) 

(54) 
(55) 
(56) 
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where  AT^  and  Ap^^  are  obtained  by  solving  simultaneously  the  algebraic 
equations  (51)  and  (52).   The  process  is  repeated  to  within  a  pre- 
determined accuracy.   The  calculations  were  repeated  for  every  degree 
of  crank  revolution  and  the  thermodynamic  properties  and  chem.ical 
composition  of  the  burnt  gases  at  exhaust  are  thus  obtained. 

5.5   The  Kinetics  of  Nitric  Oxide  Formation 

It  was  mentioned  before  that  in  the  combustion  process,  the 
formation  of  nitric  oxide  does  not  reach  the  chemical  equilibrium 
level.   In  this  section  the  reactions  and  other  factors  governing  the 
formation  of  nitric  oxide  based  on  previously  reported  work  are 
discussed,  and  the  calculation  procedure  is  outlined  in  the  next 
section. 

Combustion  processes  involving  air  as  oxidant  result  in  the 
formation  of  various  oxides  of  nitrogen:   NO,  NO,  NO   NO   etc. 
Of  all  these,  nitric  oxide  (NO)  is  predominant  in  automobile  exhausts; 
measurable  amounts  of  N0„  have  also  been  detected  in  certain  exhausts. 

The  formation  of  nitric  oxide  in  combustion  processes  is 
affected  primarily  by  the  following  three  factors: 

(a)  Tne  combustion  temperature:   Since  nitric  oxide  is  a  high-enthalpy 
species  relative  to  molecular  nitrogen  and  oxygen,  from  the  viewpoint  of 
therm.odynamics ,  the  formation  of  nitric  oxide  is  favored  by  the  existence 
of  high  combustion  temperatures. 

(b)  Air-fuel  ratio:   Tns  air-fuel  ratio  influences  the  formation  of 
nitric  oxide  through  its  effect  on  the  availability  of  excess  oxygen. 
The  formation  of  nitric  oxide  is  low  for  low  air-fuel  ratios  and 


108 


increases  to  a  maximum  for  mixtures  having  approximately  10%  excess 
air  (as  compared  to  stoichiometric  mixtures).   This  increase  is  due 
partly  to  the  increase  in  the  availability  of  oxygen  and  partly  to 
increasing  combustion  temperatures.   Further  increases  in  the  air-fuel 
ratio  decrease  the  nitric  oxide  formation  since  the  lower  combustion 
temperatures  offset  the  increased  availability  of  oxygen. 
(c)   Reaction  time  available:   This  is  a  very  important  factor  to  be 
considered  in  nitric  oxide  formation,  especially  in  transient  systems 
like  the  internal  combustion  engine. 

Early  analyses  of  the  formation  of  nitric  oxide  were   based 
on  chemical  equilibrium  analysis  assuming  that  sufficient  time  was 
available  for  the  combustion  products  to  reach  chemical  equilibrium 
[17],   If  chemical  equilibrium  is  assumed  to  exist  during  the  expansion 
stroke  also,  then  the  analysis  reveals  that  the  concentration  of  nitric 
oxide  at  exhaust  should  be  negligible  [18],  although  significant 
quantities  (several  hundred  parts  per  million  to  as  high  as  1-2  mole 
percent)  may  be  formed  at  the  completion  of  the  combustion  process. 
The  exact  amount  form.ed  depends  on  the  combustion  temperature,  the 
air-fuel  ratio  and  other  factors. 

Pvecent  investigations  have  revealed  that  although  the  assumption 
of  equilibrium  combustion  is  a  good  one  when  the  primary  combustion 
products  are  concerned,  it  does  not  portray  the  true  state  of  affairs 
as  far  as  the  formation  of  nitric  oxide  is  concerned.   Contrary  to 
that  for  the  other  species  involved,  for  nitric  oxide,  the  time 
required  to  reach  equilibrium  concentration  is  significant  compared  to 
the  time  available  in  a  typical  internal  combustion  engine  [16].   Also, 
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the  fact  that  the  rates  of  the  cheraical  reactions  involved  in  the 
formation  of  nitric  oxide  are  highly  sensitive  to  variations  in 
temperature,  coupled  with  the  drop  in  temperature  of  the  combustion 
products  following  combustion,  means  that  equilibrium  concentration 
levels  for  nitric  oxide  are  seldom,  if  ever,  reached.   The  fact  that 
the  formation  of  nitric  oxide  is  rate  controlled  and  not  governed 
by  chemical  equilibrium  considerations  has  been  demonstrated  both 
theoretically  and  experimentally  in  several  studies  involving 
different  combustion  systems  [19,  20,  21,  22,  2  3]. 

The  assumption  of  chemical  equilibrium  during  the  expansion 
stroke  is  also  a  poor  one.   The  rates  of  the  reactions  resulting  in 
the  decomposition  of  nitric  oxide  fall  to  insignificant  values  early 
in  the  expansion  stroke  due  to  the  lowering  of  temperature.   As  a 
result,  the  concentration  of  nitric  oxide  is  essentially  frozen  at 
its  value  early  in  the  expansion  stroke  [18,  19,  24,  25,  26].   It  has 
been  shov^^n  experimentally  that  the  concentration  of  nitric  oxide  at 
exhaust  is  greatly  in  excess  of  that  predicted  by  an  equilibrium 
analysis. 

The  formation  of  nitric  oxide  is  thus  controlled  by  the  reaction 
rates  during  the  high  temperature  part  of  the  process  and  the  time 
available. 

In  some  of  the  early  investigations  the  formation  of  nitric 
oxide  was,  as  before,  considered  on  an  equilibrium  analysis  cor- 
responding to  the  peak  flame  temperature  but  the  decomposition  process 
was  assumed  to  be  essentially  frozen  because  of  two  reasons:   the 
temperature  falls  rapidly  during  the  expansion  stroke  and  the  oxygen 
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atom  concentration  (involved  in  the  decomposition  reactions)  is 
relatively  low  in  the  burnt  gases  during  the  expansion  stroke. 
Such  an  analysis  would  predict  a  nitric  oxide  concentration  at 
exhaust  higher  than  actual  since,  firstly,  the  formation  process 
seJdom  reaches  the  equilibrium  level  and,  secondly,  some  decom- 
position, however  little,  does  occur  during  the  expansion  stroke. 
Newhall's  theoretical  analysis  [13]  shows  that  the  decomposition  of 
nitric  oxide  during  the  expansion  stroke  is  not  appreciable.   Other 
theoretical  analyses  [26,  27]  predict  the  same  result. 

Now  some  other  factors  affecting  the  concentration  of  nitric 
oxide  will  be  mentioned  briefly.   Ignition  timing  plays  an  important 
role  in  the  amount  of  nitric  oxide  formed,  since  it  affects  the  time 
available  for  its  formation  and  also  the  temperature.   Ihe  nitric 
oxide  emissions  increase  the  earlier  the  spark  timing  [28]. 

Due  to  increased  temperatures  at  increased  pressures,  nitric 
oxide  emissions  increase  with  increasing  pressure  (which  increases 
with  load).   Tue  same  effect  is  observed  xvith  increasing  compression 
ratios.   In  addition,  at  lower  compression  ratios,  the  fraction  of 
residual  gases  from  the  previous  cycle  is  larger  and  so  this  reduces 
the  nitric  oxide  em.issions  (same  effect  as  exhaust  gas  recirculation). 

Exhaust  gas  recirculation,  x-zater  injection,  and  increased  valve 
overlap  all  increase  the  fraction  of  inert  gases  in  the  combustion 
chamber,  thereby  reducing  flame  temperatures  and  hence  nitric  oxide 
emissions.   In  addition,  during  exhaust  gas  recirculation,  the 
additional  carbon  dioxide  and  water  present  dissociate  at  high 
temperatures,  thereby  reducing  the  peak  temperatures  and  hence  the 
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rate  of  nitric  oxide  formation,  besides  depressing  the  equilibrium 
concentration.   Water  injection  achieves  the  same  effect  through 
vaporization  and  later  dissociation.   Intake  charge  dilution  also 
reduces  the  nitric  oxide  emissions;  the  higher  the  heat  capacity  of 
the  diluent,  the  more  the  reduction  [29].   The  composition  of  the 
fuel  affects  the  nitric  oxide  emissions  due  to  the  fact  that  fuels 
with  higher  hydrogen-carbon  ratios  have  lower  peak  flame  temperatures, 
hence  lower  nitric  oxide  emissions  [17,  30]. 

It  is  thus  seen  that  the  formation  of  nitric  oxide  occurs  mostly 
in  the  burnt  gases  after  the  passage  of  the  flame.   It  has  been  shown 
[19,  21]  that  the  formation  of  nitric  oxide  in  post  flam.e  gases  may  be 
accurately  predicted  by  the  following  simple  mechanism  (known  as  the 
Zeldovich  mechanism) : 


N  +  0  =  NO  +  N  [1] 

0  +  N  =  NO  +  0  [2] 


"£""  ^  kif[N2][0]  -  k^^[NO][N]  +  k2^[02][N]  -  k2^[N0][0]   (57 


) 


Since    the    rate    of  nitric   oxide    formation  is   slow   compared   to   the   rates 
of    the   other   chemical   reactions,    the   equilibrium  concentrations   of  N, 
N        0,    0^    may   be    used    in    the    above    equation.      Using   this    theoretical 
procedure,    excellent    agreement  with   the   experimentally  measured 
concentrations    of  nitric   oxide    in  hydrogen-air   and  propane-air    [19  ,   21 ] 
mixtures  has   been  obtained    (Figs.    5.5,    5.6). 
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Experimental   and   theoretical  nitric   oxide 
formation   rates    for   combustion  of  hydrogen- 
air  mixtures    (Ref,    19) 
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formation  rates  for  combustion  of  propane- 
air  mixtures  (Ref.  21) 
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Besides  the  above  two  reactions ,  some  other  reactions  involved 
in  the  formation  of  nitric  oxide  are  3,  4,  6,  7,  and  8  of  Table  5.1  and 

N  +  OH  =  NO  +  H  [18] 

NO  +  0   =  NO  +  NO  [19] 

NO  +  NO   =  NO  +  NO  +  0  [20] 

Although  some  of  these  reactions  have  large  rate  constants ,  they  in- 
volve species  present  in  very  small  concentrations  and  hence  their 
contribution  towards  the  formation  is  negligible  compared  to  the 
reactions  (1)  and  (2).   Under  fuel- rich  conditions  reaction  (18)  may 
be  significant;  however,  there  is  some  controversy  as  to  the  reported 
value  of  the  rate  constant  for  this  reaction.   Under  fuel-lean  conditions 
reaction  (3)  may  have  a  significant  contribution. 

5.6   Details  of  Calculations  for  Nitric  Oxide  Formation 

Equation  (57)  expressing  the  rate  of  formation  of  nitric  oxide 
was  integrated  for  a  time  interval  corresponding  to  one  degree  of  crank 
revolution,  and  the  concentration  of  nitric  oxide  thus  obtained  was  used 
in  the  calculations  involving  the  flame  element,  while  the  concentrations 
of  the  other  species  were  obtained  from  a  chemical  equilibrium  analysis 
(Section  5.4.3).   For  the  burnt  gases,  the  nitric  oxide  form-aticn 
equation  was  again  integrated  over  one  degree  of  crank  revolution  at 
each  step  of  the  computational  process. 

Thus  the  concentration  of  nitric  oxide  at  the  end  of  the 
combustion  process  was  obtained  on  the  basis  of  a  rate-dependent 
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analysis.   The  final  concentration  of  nitric  oxide  at  any  point 
depends  on  the  temperature  at  that  point  and  the  time  interval 
between  the  arrival  of  the  flame  at  that  point  and  the  time  at 
which  temperatures  in  the  expansion  stroke  fall  to  values  below  which 
the  rate  of  formation  is  negligible.   Since  the  formation  of  nitric 
oxide  begins  earlier  at  locations  nearer  the  point  of  ignition,  larger 
amounts  of  nitric  oxide  are  formed  at  such  locations.   This  is  further 
augmented  by  the  fact  that  temperatures  are  also  higher  at  locations 
nearer  the  point  of  ignition  and  the  fact  that  the  rates  of  the 
reactions  in  question  are  sensitive  to  temperature.   This  has  been 
confirmed  experimentally  [ 20] . 

Thus,  in  reality,  at  the  end  of  the  combustion  process  there 
is  a  gradient  of  temperature  and  nitric  oxide  concentration  in  the 
combustion  chamber  with  higher  temperatures  and  higher  concentrations 
nearer  the  point  of  ignition,  which  in  som.e  cases  approach  equilibrium 
values  [31j.   In  this  analysis  transport  effects  were  considered  to 
be  absent.   Even  if  considered,  the  transport  processes  of  diffusion 
and  thermal  conduction  (or  convection  and  radiation)  are  unable  to 
equalize  the  concentration  and  temperature  gradients.   In  the  present 
work,  a  bulk  analysis  has  been  made  wherein  the  burnt  gases  in  the 
combustion  chamber  are  assumed  to  be  at  one  uniform  'pseudotemperature ' 
and  'pseudoconcentration'  obtained  by  an  averaging  process. 

5.7   Injection  of  Water-Alcohol  Mixtures 

Now  the  modifications  to  be  made  in  the  cycle  analysis  due  to 
the  injection  of  water-alcohol  mixtures  will  be  considered.   Mien  this 
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mixture  is  injected  into  the  intake  manifold,  the  analysis  is  the 
same  as  before,  except  that  the  concentrations  of  the  species  in  the 
incoming  fuel-air-alcohol  mixtures  have  to  be  suitably  adjusted.   IVhen 
the  alcohol-water  mixture  is  injected  in  the  middle  of  the  combustion 
process,  certain  modifications  have  to  be  made.   These  v/ere  considered 
on  the  basis  of  three  different  models  which  are  discussed  below.   All 
the  models  involve  approximations  in  varying  degrees.   The  ignorance 
of  the  rate  dependent  processes  (chemical  kinetics)  occurring  after 
injection  prevents  a  thorough  realistic  analysis  from  being  made. 
In  the  absence  of  experimental  data  the  accuracy  of  the  different 
models  could  not  be  compared  or  checked. 

In  all  the  three  cases  it  is  assumed  that  injection  occurs 
at  the  beginning  of  a  combustion  step  (Section  5.4.3)  when  the 
engine  cylinder  contains  only  the  burnt  and  unburnt  fractions.   The 
changes  occurring  due  to  injection  are  assumed  to  be  completed  before 
the  combustion  step  is  initiated.   In  practice,  the  processes  would 
occur  at  finite  rates  but  since  reaction  times  for  hydrocarbons  at 
high  temperatures  are  small,  the  reactions  m.ay  be  assumed  to  be 
instantaneous  to  a  first  approximation.   Tlie  injected  mixture  is 
assumed  to  be  divided  between  the  burnt  and  unburnt  fractions  in  the 
ratio  of  their  respective  volum.es. 

On  the  unburnt  side,  it  is  assumed  that  no  chemical  reactions 
occur,  and  that  the  result  is  a  homogeneous  mixture  at  a  uniform 
pressure  and  tem.perature .   The  concentrations  of  the  species  in  this 
mixture  were  found  from  a  relation  similar  to  equation  (2)  and  the 
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temperature  was  obtained  from  an  iterative  process  similar  to  that 
described  by  equations  (36-40). 

The  three  models  for  considering  the  sequence  of  events  on 
the  burnt  side  are: 

(1)  Constant  pressure  equilibrium:   This  assumes  that  the  injected 
mixture  and  the  burnt  gases  are  brought  to  a  chemical  equilibrium  at 
constant  pressure.   The  equilibrium  composition  was  obtained  by  the 
method  outlined  in  Appendix  I  and  the  temperature  was  obtained  by  an 
iterative  process  similar  to  that  described  by  equations  (3-8)  utilizing 
the  condition  that  the  total  enthalpy  remains  constant  (constant  pressure 
adiabatic  process). 

(2)  Constant  volume  equilibrium:   This  is  the  same  as  above  except 
that  now  equilibrium  is  achieved  at  constant  volums.   Tlie  iteration 
for  the  temperature  was  based  on  the  condition  that  the  total  internal 
energy  (instead  of  the  total  enthalpy)  remains  constant,  and  equations 
(36-40)  xjere  used. 

(3)  Equilibrium  decomposition  of  droplets:   This  assum.es  that  the 
injected  mixture  is  present  in  the  burnt  gases  in  the  form  of  uniform 
droplets  and  decomposes  in  the  presence  of  the  ox>-gen  available  to  a 
chemical  equilibrium  composition.   Constant  pressure  decomposition  is 
assumed  and  so  the  temperature  and  equilibrium  composition  of  the 
resulting  products  were  obtained  as  in  the  first  stage  of  the  combustion 
step.   These  resulting  products  were  then  mixed  with  the  burnt  gases 
(assuming  no  further  chemical  reactions)  and  the  resulting  temperature 
found  on  the  basis  of  a  constant  pressure  or  constant  volum.e  assumption 
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as  in  the  mixing  step  described  before;  equations  (36-40)  or  (41-45) 
were  used  respectively. 

For  all  the  three  cases  the  final  step  was  to  restore  the 
total  system  to  its  volume  before  the  injection  process,  and  this 
was  done  using  the  set  of  equations  (33-35). 
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CHAPTER  VI 
RESULTS  AND  DISCUSSIONS 

6.1   The  Coiribustion  Process 

In  this  section  the  results  obtained  from  the  inatheiTiatica.l 
analysis  of  the  combustion  process  as  described  in  Chapter  IV  are 
presented. 

6.1.1   Concentration  and  temperature  profiles 

The  solution  of  tlie  set  of  partial  differential  equations 
governing  one-dimensional  unsteady  flame  propagation  yields  the 
distribution  of  the  species  concentrations  and  temperature  through  the 
flame  zone.   The  results  are  presented  graphically  in  the  following 
pages. 

Fig.  6.1  shows  the  temperature  profile  at  successive  instants 
of  time  and  shows  the  development  of  the  steady  state  profile  from 
the  starting  one.   It  is  clear  that  by  choosing  an  initial  profile 
closely  resembling  the  steady  state  one,  the  computation  can  be  cut 
down  considerably. 

Figs.  6.2,  6.3,  6.4,  and  6.5  give  the  distribution  of  the  species 
mole  fractions  through  the  flame  zone.   It  is  observed  that  the  concentra- 
tions of  the  four  major  species  C0„ ,  H_0  (major  products),  0„ ,  and  CH^ 
(major  reactants)  either  increase  or  decrease  continuously  from  the 
cold  to  the  hot  boundary. 
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the  stead}'  state  profile 
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Fig.  6.3.   Concentration  profiles  for  the  species  H„ , 
OH,  NO,  NO^,  and  N^O 
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Fig.  6.5.   Concentration  profiles  for  the  species  CH^, 
CHO,  and  N 
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The  methane  and  oj:ygen  profiles  show  a  continuous  decrease 
from  the  cold  to  the  hot  boundary  showing  that  they  are  being 
continuously  used  up.   The  carbon  dioxide  and  vjater  profiles  exhibit 
an  opposite  trend.   The  disappearance  of  methane  is  due  to  the  chain 
propagation  reactions: 

1.  CH,  +  OH  =  CH^  +  H^O 

4  3     2 

2.  CH,  +  H  =  CH„  +  H„ 

4  5  1 

3.  CH,  +  0  =  CH^  +  OH 

4  J 

A  study  of  the  concentration  profiles  reveals  that  at  the 
point  where  the  methane  disappearance  starts,  the  concentration  of 
the  K-atoms  is  very  low  compared  to  the  concentrations  of  the 
OK-radicals  and  the  0-atoms,  so  that  the  second  reaction  may  be 
neglected  here.   This  reaction  may,  however,  be  important  for  fuel-rich 
flam.es.   Reaction  3  is  slower  than  reactions  1  and  2  [1]  and  so  reac- 
tion 1  is  the  dominant  one  in  the  disappearance  of  methane  in  the 
initial  stages.   Tnis  reaction  also  results  in  the  formation  of  water 
and  generates  CH„-radicals.   In  the  reaction  mechanism  given  in 
Table  4.1,  all  these  reactions  have  been  considered,  since  reaction  3, 
and  to  a  lesser  extent  reaction  2,  does  have  some  contribution  in  the 
later  stages  of  the  process.   Another  reaction  considered  in  the  dis- 
appearance of  methane  is: 

4.  CH,  +  M  =  CH^  +  H  +  M 

4  J 

This  reaction  has  a  significant  contribution  in  some  regions  of  the 
flame.  However,  when  the  concentration  of  the  radicals  has  reached 
a  high  value  in  the  primary  reaction  zone,  the  reverse  reaction 
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predominates  the  forward  one,  resulting  In  the  recombination  of  the 
CH„-radlcals  and  the  H-atoms.   This  being  a  three-body  reaction,  Its 
rate  is  small. 

The  generation  of  OH-radicals  required  for  reaction  1  can  occur 
in  the  following  ways: 

5.  H  +  0   =  OH  +  0 

6.  0  +  H   =  OH  +  H 

7.  0  +  HO  =  OH  +  OH 

8.  H  +  HO  =  H  +  OH 

A  study  of  the  relative  concentrations  of  the  various  species 
Involved  reveals  that  reaction  5  is  mainly  responsible  for  the  genera- 
tion of  OH-radlcals  and  this  reaction  also  results  in  the  disappearance 
of  oxygen.   Reaction  6  is  of  secondary  importance;  however,  it  could 
play  a  major  role  In  fuel-rich  flames. 

The  production  of  carbon  monoxide  in  the  initial  stages  of  the 
flame  occurs  due  to  the  reaction; 

9.  OH  +  CHO  =  CO  +  HO 

where  the  CHO-radicals  are  obtained  from  the  reactions  involving  the 
CH  -radicals  and  the  oxygen  atoms  or  molecules: 

10.  CH  +  0^  =  H2O  +  CHO 

11.  CH^  +  0   =  H^  +  CHO 

Finally,  some  recombination  reactions  which  play  a  role  in  the 
latter  stages  of  the  process  are: 

12.  CO  +  H  +  M  =  CHO  +  M 

13.  H  +  OH  +  M  =  HO  +  M 

14.  H  +  0   +  M  =  OH  +  M 
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The  concentration  profiles  for  the  intermediate  species  exhibit 
a  behavior  which  is  in  general  different  from  that  for  the  main  species. 
The  free  radicals  have  a  maximum  concentration  somewhere  within  the 
flame  zone,  the  exact  location  being  a  complicated  function  of  the 
diffusion  and  the  chemical  rate  of  production  of  species.   The  concentra- 
tions of  the  free  radicals  in  the  flame  zone  are  considerably  higher 
(quite  often  by  a  few  orders  of  magnitude)  than  their  equilibrium 
concentrations  at  the  flame  boundaries.   The  free  radicals  play  an 
important  role  in  the  flame  reactions — chain  initiation  and  chain 
branching,  and  this  accounts  for  the  rapidity  of  the  reactions. 

In  the  primary  reaction  zone,  the  bimolecular  reaction  steps 
1-8  are  of  importance  and  are  responsible  for  the  rapid  generation  of 
large  quantities  of  the  free  radicals.   These  are  the  chain  initiation 
and  chain  branching  steps.   The  removal  of  radicals  occurs  by  chain 
breaking  steps.   These  are  mainly  the  three-body  reactions  12-14.   ITie 
relative  slowness  of  the  chain  breaking  steps  accounts  for  the  radical 
concentrations,  in  excess  of  the  chemical  equilibrium  concentrations 
being  generated  in  the  initial  stages  of  the  process  before  the  three- 
body  reactions  become  effective  in  reducing  the  radical  concentrations. 

Experimental  data  on  the  concentration  and  temperature  profiles 
have  been  reported  for  low  pressure  (0.1  atm. )  methane-air  flames  [2]. 
Tne  concentrations  of  the  major  species  CH, ,  0    ,    CO   CO,  and  K„0  have 
been  experimentally  measured  and  these  show  good  agreement  with  the 
theoretical  results  obtained  here. 
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6.1.2   Unity  Lewis  number  approximations 

When  the  Lewis  number  pC  D . . /A  is  equal  to  unity  for  all  the 

P  ij 

species  pairs  in  the  mixture  of  mean  specific  heat  C   and  thermal 
conductivity  A,  the  diffusion  and  conduction  fluxes  cancel  everywhere 
in  the  flame  zone  and  hence  the  enthalpy  is  constant  through  the  flame 
zone.   Klein  [3]  has  shown  that  for  a  complex  flame  if  the  heat  release 
can  be  related  to  the  consumption  of  some  single  'fuel',  then  the  Lewis 
number  of  this  fuel  and  the  corresponding  oxidizer  will  govern  the 
behavior  and  the  Lewis  numbers  of  the  other  species  pairs  would  have  a 
minor  effect  on  the  overall  behavior.   Under  these  conditions  the  actual 
complex  flame  reduces  to  a  simple  one  with  A->B  type  of  behavior,  so  that 
the  concentration  of  the  fuel  is  a  linear  function  of  the  temperature: 


T  -  T 


fuel    '  fuel  o  T  -  T 
f    c 


where   T   =  initial  temnerature  corresponding  to  (X-   ^ ) 
o  to      fuel  o 

T  =  final  temperature  at  which  all  the  fuel  is  consumed 

The  methane-air  flame  has  two  reaction  zones,  one  in  which 
methane  is  burnt  and  the  other  in  which  carbon  m.onoxide  is  burnt  so 
that  these  two  may  be  considered  as  the  fuel  in  their  respective  zones. 
Tne   T   for  carbon  monoxide  and  methane  are  approximately  1960  degrees  K 
and  630  degrees  K  respectively  (Fig.  6.1),  while  the  T  values  are 
2390  degrees  K  and  1950  degrees  K  respectively. 

A  plot  of  the  concentration  versus  the  dimensionless  tem.perature 
ratio  6  =  (T  -  T  )/(T  -  T  )  shows  approximately  linear  behavior  which 


132 


shows  that  the  assumption  of  unity  Lewis  number  may  be  used  to  obtain 
an  approximate  solution  for  the  flame  under  consideration  (Fig.  6.6), 

The  experimental  results  obtained  by  Fristrom  and  Westenberg  [4] 
for  a  0.1  atm.  methane-oxygen  flame  are  shown  in  Fig.  6.7  for  the  sake 
of  comparison.   The  agreement  with  the  present  results  is  excellent. 

Fig.  6.8  shows  the  steady  state  temperature  profiles  for  the 
cases  of  unity  and  non-unity  Lewis  numbers.   The  closeness  of  the  Lewis 
number  to  unity  is  clearly  borne  out.   The  actual  values  of  the  Lewis 
numbers  of  methane-oxygen  and  carbon  monoxide-oxygen  are  1.25  and  1.11 
respectively,  and  so  the  behavior  noted  above  is  not  surprising. 

6.1.3  Two  stage  hydrocarbon  combustion 

Experimental  data  on  propane-air  flames  [5,  6]  and  other  studies 
[7]  have  revealed  that  the  combustion  of  hydrocarbons  occurs  in  two 
stages.   In  the  first  stage  the  hydrocarbon  is  converted  to  carbon 
monoxide  and  water,  and  in  the  second  stage  carbon  monoxide  is  burnt 
to  carbon  dioxide.   The  results  obtained  here  corroborate  this  fact. 
Firstly,  the  unity  Lewis  number  approximation  discussed  above  points 
to  the  fact  that  combustion  might  be  occurring  in  two  different  stages. 
A  look  at  the  concentration  profiles  tends  to  confirm  this. 

In  the  initial  stages  of  the  flame,  methane  is  consumed  and 
carbon  monoxide  and  water  are  produced.   At  the  point  where  the  methane 
concentration  has  been  reduced  to  negligible  values,  the  carbon  monoxide 
concentration  attains  its  maximum  va].ue  and  the  x^/ater  concentration  is 
close  to  its  final  value  at  the  hot  boundary.   After  this  the  carbon 
monoxide  concentration  decreases  witVi  a  corresponding  increase  in  the 


133 


0.10 


0.08 


^  0.06 

o 

•H 

o 

u 

!°0.04 


0.02 


0.00 


1.0 


e  = 


(T-T^)/(T^- 


^o> 


Fig.  6.6.   Test  of  unity  Lewis  number  approximation  for  methane 
and  carbon  monoxide  in  a  methane-air  flame.   Linear 
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carbon  dioxide  concentration,  and  this  is  the  portion  of  the  flame 
where  most  of  the  carbon  dioxide  is  produced.   The  disappearance  of 
carbon  monoxide  is  mainly  due  to  the  reaction 

CO  +  OH  =  CO  +  H 

This    confirms   the   experimentally   observed  phenomenon   of   two-stage 
hydrocarbon   combustion. 

6.1.4   Closing  remarks 

Some  results  obtained  for  the  methane-air  flame  have  been 
discussed  in  the  preceding  sections.   Tlie  main  objective  of  this  study 
has  been  to  present  a  general  computational  scheme  to  obtain  a  detailed 
description  of  the  flame  zone.   This  has  been  done  and  the  results 
obtained  here  have  been  compared  with  the  limited  experimental  data 
available  and  the  validity  of  the  computational  method  has  been 
demonstrated.   No  attempt  has  been  made  to  draw  conclusions  for 
hydrocarbon  flames  in  general,  since  the  results  obtained  would  not 
warrant  such  an  attempt. 

As  experim.ental  data  on  other  flames  are  gathered,  the  method 
may  be  used  to  gain  more  insight  into  the  various  phenomena  occurring 
in  the  flame  zone  and  obtain  a  better  understanding  of  the  chemical 
kinetics  of  the  flame  reactions.   The  latter  is  perhaps  the  most 
important  use  of  flame  structure  data. 

The  improvements  made  in  this  study  over  previous  studies  m.ay 
be  summed  up  as: 
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(1)     A  more  accurate  evaluation  of  the  transport  and  thermodynamic 
properties,  the  evaluation  being  limited  by  the  lack  or  scarcity  of 
experimental  data  for  the  transport  coefficients  ,  especially  at  high 
temperatures. 

(ii)    A  realistic  consideration  of  the  chemical  processes.   The  concept 
of  a  single  overall  reaction  governing  the  chemical  processes  has  been 
discarded  as  unreal.   Instead,  a  set  of  elementary  reaction  steps  has 
been  considered  to  govern  the  production  or  destruction  of  the  chemical 
species,  and  thereby  the  important  role  played  by  free  radicals  in  the 
flame  reactions  has  been  taken  into  account. 

6.2  Mathematical  I-lodel  for  the  Otto  Cycle 

The  mathematical  model  described  in  Chapter  V  was  used  to  obtain 
the  state  properties  and  the  species  concentrations  at  every  stage  in 
the  cycle.   The  engine  specd.fications  corresponding  to  which  calculations 
were  carried  out  are  given  in  Appendix  III  and  represent  a  standard  144 
CID  six-cylinder  engine.   The  choice  of  engine  size  and  configuration  was 
completely  arbitrary.   A  specific  choice  was  necessary  to  demonstrate  the 
procedure  and  make  comparisons . 

A  graphical  representation  of  the  results  of  interest  is  given 
in  the  following  pages,  v/hile  a  typical  computer  output  is  provided  in 
Appendix  III.   These  results  are  for  the  case  of  a  stoichiometric  mixture 
of  methane  and  air. 

6.2.1   The  combustion  process 

The  combustion  process  is  initiated  at  30  degrees  before  top 
dead  center  and  completed  10  degrees  after  top  dead  center.   Fig.  6.9 
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shows  the  rise  of  pressure  and  the  temperatures  of  the  burnt  and  the 
unburnt  fractions  during  this  process.   It  is  observed  that  all  these 
three  quantities  increase  rapidly  during  the  initial  stages  of  the 
combustion  process.   The  increase  is  associated  mainly  with  the  energy 
released  by  the  burning  of  the  fuel  and  to  a  smaller  extent  due  to  the 
compression  of  the  gases  to  a  smaller  volume  until  top  dead  center  is 
reached.   At  Lop  dead  center,  most  of  the  fuel  has  undergone  combustion 
(more  than  98%).   Tne  energy  released  by  the  burning  of  the  remainder 
of  the  fuel  is  not  enough  to  prevent  the  drop  in  the  burnt  and  unburnt 
gas  temperatures  and  the  pressure  after  top  dead  center.   This  drop  is 
partly  due  to  the  increasing  volume  as  the  expansion  stroke  comm.ences 
and  partly  due  to  the  heat  transfer  to  the  relatively  cool  combustion 
cham.ber  walls.   The  burnt  and  the  unburnt  gas  temperatures  and  the 
pressure  all  reach  their  peak  values  at  top  dead  center  after  which 
they  start  dropping  slowly. 

Fig.  6.10  shows  the  concentration-time  history  of  nitric  oxide 
in  the  burnt  gases  during  the  combustion  process.   In  the  initial  stages 
the  rate  of  formation  of  nitric  oxide  is  slow  because  of  the  relatively 
low  burnt  gas  temperature  and  the  correspondingly  low  concentration  of 
ox>'gen  and  nitrogen  atoms.   However,  before  combustion  has  proceeded 
for  long,  the  rate  of  formation  and  the  concentration  of  nitric  oxide 
increase  rapidly  due  to  rapidly  increasing  burnt  gas  temperatures  and 
radical  concentrations.   Nitric  oxide  formation  continues  even  after 
top  dead  center  is  reached  but  the  process  is  considerably  slower, 
partly  because  the  burnt  gas  temperature  starts  falling,  and  partly 
because  the  concentration  has  reached  a  good  fraction  of  the  chemical 


139 


1100 


1000  - 


900 


800  - 


700  - 


600 


-30 


Fig.  6.9 


-20         -10  0  10 

CRANK  AI^'GLE  (DEGREES  FROM  TDC) 

Pressure  and  temperature  histories  during 

combustion  process 

T  :   Burnt  gas  temperature 


T  :   Unburnt  gas  temperature 

U  or 

p:   Pressure  (kg. /sq. cm. ) 


J   (degrees  K) 


140 


-30 


10 


-20         -10  0 

Crank  angle  (degrees  from  TDC) 
Fig.  6,10.   Nitric  oxide  formation  during  conibustion 


141 


equilibrium  concentration.   The  nitric  oxide  concentration  for  this 
case  was  found  to  be  a  maximum  of  slightly  over  1  mole  percent  at 
about  9  degrees  after  top  dead  center;  the  corresponding  burnt  gas 
temperature  and  pressure  V7ere  2918  degrees  K  and  74,70  kg./sq.cm.  , 
respectively,  while  the  peak  values  of  the  tem.perature  and  pressure 
were  respectively  2963  degrees  K  and  80  kg./sq.cm.,  both  at  top  dead 
center. 

The  mathematical  model  also  yields  the  concentration  time 
histories  of  the  other  species  for  the  duration  of  the  combustion 
process  and  the  changes  in  the  state  properties:   ratio  of  specific 

heats  (C  /C  ),  the  soecific  heat  at  constant  pressure,  enthalpy, 
p   V  ' 

internal  energy,  entropy,  and  the  Gibbs  free  energy.   Tliese  may  be 
found  in  the  computer  output  in  Appendix  III. 

6.2.2  Tae   expansion  process 

The  computed  pressure-time  and  temperature-time  histories 
for  the  expansion  process  are  presented  in  Fig.  6.11.   It  is  observed 
that  both  the  pressure  and  the  temperature  decrease  rapidly,  the 
former  more  rapidly  than  the  latter,  during  the  initial  stages  of  the 
expansion  process.   The  decrease  is  mainly  due  to  the  increasing 
volume  and  to  a  smaller  extent  due  to  the  heat  transfer  from  the  hot 
gases  and  the  fact  that  radical  recombination  reactions  decrease  the 
effective  total  number  of  moles  of  gases  present.   The  recombination 
reactions  also  give  off  energy,  partly  counteracting  the  effects  of 
the  heat  losses  to  the  surroundings.   The  effect  of  the  chemical 
reactions  on  the  pressure  and  temperature  is  of  minor  importance 


142 


( 'luo  "bs/ "S^i)    aanssaJ^j 


o 

OD 


o 


o 


o 


0) 

^ 

4J 

^ 

o 

U-l 

o 

CO 

w 

<-i 

4-1 
O 

p. 

,_^ 

0) 

u 

e 

e 

•H 
4-) 
1 

!-i 

1 

0) 

^1 

J-) 

3 

o 

in 

4J 

ro 

tu 

CO 

r-t 

U 

cn 

<u 

Q) 

(i 

Q) 

B 

!-i 

0) 

to 

4J 

(D 

T) 

XI 

s_^ 

S 

G 

CS 

y 

0) 

O 

o 

M 

0) 

^ 

cr\ 

bO 

B 

4J 

a 

•H 

cf; 

C3 

4-1 

1 

C 

^ 

Q) 

O 

a 

^1 

•H 

C3 

3 

tfi 

J-J 

tn 

*-< 

u 

cn 

rt 

QJ 

Cu 

^^ 

^ 

O 

Cl, 

Q) 

in 
o 

• 
rH 

• 

60 
•H 

<-t 

fn 

o 
o 
o 


O 
O 
vD 


o 
o 

CN 
CM 


O 
O 
tX3 


O 

o 


()I    S33J:§3p)    aan^Basduisx 


143 


because  of  either  one  or  both  of  the  following  two  reasons.   Firstly, 
the  net  energy  release  or  absorption  due  to  the  chemical  reactions 
is  small  compared  to  the  energy  associated  with  the  work  done  during 
the  expansion  stroke.   This  is  due  to  the  comparatively  low  heats 
of  reaction.   Secondly,  the  rates  of  the  reactions  are  slow  especially 
as  the  temperature  drops  to  relatively  low  values  in  the  later  parts 
of  the  expansion  stroke.   The  pressure  curve  almost  levels  off  around 
130  degrees  after  top  dead  center,  while  the  temperature  continues  to 
drop,  though  at  an  increasingly  slower  rate,  up  to  the  end  of  the 
expansion  stroke, 

Tlie  concentration-time  histories  of  the  various  species  are 
presented  in  Figs.  6.12,  6.13,  6.14,  6.15,  6.16,  and  are  discussed 
below. 

Nitric  oxide  concentrations  freeze  out  early  during  the 
expansion  stroke,  agreeing  with  resiilts  observed  in  previously 
reported  work  [8,  9].   The  reason  for  this  is  that  the  rates  of  the 
reactions  involved  in  the  decomposition  of  nitric  oxide  are  extremely 
sensitive  to  temperature  and  fall  to  negligible  values  as  the  temper- 
ature rapidl]/  falls  during  the  expansion  stroke.   In  this  case  the 
nitric  oxide  concentration  V7as  found  to  freeze  out  when  the  temperature 
had  fallen  to  about  2150  degrees  K  at  about  50  degrees  after  TDC. 
Other  reasons  contributing  to  this  phenomenon  are  the  decrease  in  the 
concentrations  of  nitrogen  and  ox^'gen  atoms  which  are  involved  in  the 
decomposition  reactions.   As  a  result,  the  nitric  oxide  concentration 
at  exliaust  is  considerably  higher  than  that  if  chemical  equilibrium 
were  to  prevail. 
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Fig.  6.12.   Concentration-time  histories  of  the  species  CO, 
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Fig.  6.13.   Concentration-time  historj'  of  species  N  during 
the  expansion  stroke 
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Fig.  6.14.   Concentration-time  history  of  species  0   during 
the  expansion  stroke 
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5.15.   Concentration-tiEB  histories  of  the  species  H„ 
and  NO   during  the  expansion  stroke 
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A  comparison  of  the  concentrations  of  the  other  species  at 
exhaust  with  the  corresponding  chemical  equilibrium  values  reveals 
that  the  concentration-time  histories  are  again  dependent  on  the  finite 
rate  of  the  chemical  reactions  and  do  not  follow  the  chemical  equilibrium 

pattern. 

Carbon  monoxide  is  formed  in  significant  amounts  (mole  fraction 
about  0.01)  during  the  coiAustion  process.   Much  of  this  gets  converted 
to  carbon  dioxide  during  the  expansion  process  as  shown  by  Fig.  6.12; 
the  reaction  CO  +  OH  =  H  +  CO^  is  mainly  responsible  for  the  oxidation 
of  carbon  monoxide.   During  the  initial  stages  the  conversion  rate  is 
fast  enough  that  the  concentration  of  carbon  monoxide  closely  follows 
its  chemical  equilibrium  value.   During  the  latter  stages,  however, 
due  to  the  reduction  in  the  concentration  of  the  radical  OH  and  the 
drop  in  temperature,  the  carbon  monoxide  concentration  deviates 
more  and  more  from  its  chemical  equilibrium  value.   Were  chemical 
equilibrium  to  prevail  at  exhaust,  the  carbon  monoxide  concentration 
would  be  much  lower  than  its  actual  value. 

The  concentration  of  all  the  radicals  (H,  N,  0,  OH)  decrease 
during  the  expansion  process  as  shown  in  Figs.  6.12  and  6.13.   The 
processes  are  again  rate-limited,  resulting  in  the  concentrations  of 
these  radicals  being  higher  than  those  corresponding  to  chemical 
equilibrium  values, 

Tne  concentrations  of  carbon  dioxide  and  water  increase 
rapidly  in  the  initial  stages  (Fig.  6.16);  carbon  dioxide,  because 
of  the  oxidation  of  carbon  m.onoxide ,  and  water,  because  of  radical 
recombination.   The  two,  however,  level  off  during  the  latter  stages 
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of  the  expansion  process  (approximately  about  90  degrees  after  TDC) 
due  to  the  fall  in  the  rates  of  the  reactions  leading  to  their 
production. 

The  variations  of  the  following  thermodynamic  properties 
during  the  expansion  stroke  may  be  found  in  the  computer  output 
provided  in  Appendix  III:   the  specific  heat  at  constant  pressure, 
the  ratio  of  specific  heats  (C  /C  ),  enthalpy,  entropy,  internal 
energy,  and  the  Gibbs  free  energy. 

Finally,  the  results  for  the  use  of  a  secondary  mixture  of 
water  and  alcohols  are  presented.   The  choice  of  this  particular 
secondary  mixture  was  arbitrary,  and  the  model  can  be  used  to  study 
the  use  of  different  substances  as  secondary  mixtures,  as  long  as 
the  assumptions  made  in  the  analysis  (Sec.  5.3)  are  not  violated 
and  provided  thermodynamic  data  for  the  substance (s)  in  question 
are  available. 

6.2.3   Water-alcohol  injection 

Firstly,  the  effect  of  injecting  the  mixture  into  the  intake 
manifold  is  considered.   Tne  results  are  summarized  In  Table  6.1. 
In  all  cases  ,  the  composition  of  the  incom.ing  fuel-air  mixture  was 
adjusted  so  that  it  was  stoichiometric  with  respect  to  (fuel+ 
methanol+ethanol) . 

It  is  seen  that  all  cases  result  in  a  reduction  of  nitric 
oxide  and  carbon  monoxide  concentrations  at  the  point  where  combustion 
is  completed,  showing  that  the  process  of  formation  of  these  pollutants 
is  definitely  retarded.   However,  the  concentrations  of  the  pollutants 
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at  exhaust  are  of  more  practical  interest.   Here  it  is  observed  that 
the  higher  the  concentration  of  carbon  monoxide  after  combustion  is 
completed,  the  lower  is  its  concentration  at  exhaust.   In  other  words, 
although  lower  amounts  of  carbon  monoxide  are  formed  during  the 
combustion  process  due  to  the  induction  of  the  mixture  of  water  and 
alcohols,  the  conversion  of  carbon  monoxide  to  carbon  dioxide  during 
the  expansion  stroke  is  also  correspondingly  slowed  dox^m.   The  result 
is  that  the  final  concentration  of  carbon  monoxide  at  ejdiaust  is  more 
than  that  for  the  case  of  no  induction.   Let  us  look  at  the  reasons 
for  this.   Firstly,  the  induction  process  causes  a  general  lowering 
of  temperature  throughout  the  cycle.   This  affects  the  rates  of  the 
reactions  during  the  expansion  stroke,  reducing  the  rate  of  conversion 
of  carbon  monoxide  to  carbon  dioxide.   Secondly,  due  to  the  lowered 
temperatures,  the  concentration  of  the  hydroxyl  radical  after 
combustion  is  completed  is  lower  and  this  adversely  affects  the 
conversion  of  carbon  monoxide  to  carbon  dioxide,  which  is  caused 
mainly  by  the  reaction  CO  +  OH  =  CO  +  H. 

For  the  case  of  nitric  oxide,  the  formation  during  the 
combustion  is  definitely  retarded  due  to  the  induction.   The  main 
reason  for  this  is  the  lowered  peak  temperatures.   However,  due  to 
the  lowered  temperatures ,  the  nitric  oxide  concentration  freezes  out 
earlier  during  the  expansion  stroke,  so  that  the  reduction  of  nitric 
oxide  at  exhaust  is  not  proportional  to  the  reduction  in  the  amount 
in  which  it  is  formed  during  combustion. 

Now  a  comparison  of  the  relative  efficiency  of  the  three 
injected  substances,  methanol,  ethanol,  and  water,  is  made.   It  is 
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found  that  the  injection  of  water  alone  is  most  effective  in  reducing 
the  nitric  oxide  emissions.   The  more  the  amount  of  water  injected  the 
more  the  reduction.   However,  the  carbon  monoxide  emissions  increase 
at  the  same  time.   Methanol  and  ethanol  have  about  the  same  effect; 
ethanol  causes  a  slight  increase  in  the  nitric  oxide  emissions  and 
a  slight  decrease  in  the  carbon  monoxide  emissions  as  compared  to 
methanol.   The  difference  is  not  significant;  less  than  5%  even  at 
100%  induction. 

Now  the  effect  of  injecting  the  mixture  in  the  middle  of  the 
combustion  process  is  considered.   Three  different  models  for  the 
analysis  of  this  process  were  presented  in  Chapter  V.   llie  results 
obtained  using  models  (1)  and  (2)  were  not  significantly  different 
(less  than  0.05%  difference  for  the  species  concentrations  and  the 
state  properties).   Tne  results  of  model  (3)  were,  however,  quite 
different  from  those  of  the  previous  two.   In  the  absence  of  experi- 
mental data,  it  is  not  possible  to  verify  the  accuracy  of  these 
calculations . 

K'hen  this  mixture  is  injected  in  the  middle  of  the  combustion 
process,  the  air-fuel  ratio  can  be  adjusted  so  as  to  be  overall 
scoicUiometric  at  the  end  of  the  injection  process.   For  this  case  the 
mixture  would  be  slightly  fuel- lean  until  injection  occurred,  and  this 
is  approximately  in  the  neighborhood  of  the  air-fuel  ratio  v/hich  causes 
maxim.um  nitric  oxide  form.ation.   Thus  the  result  is  the  formation  of 
nitric  oxide  in  large  amounts  until  injection  occurs.   Even  after 
injection  it  was  found  that  the  alcohols  were  not  very  effective  in 
reducing  the  formation  of  nitric  oxide.   Water  was  effective  since  it 
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caused  lowered  peak  temperatures,  but  a  side  effect  was  largely 
increased  carbon  monoxide  emissions. 

The  other  alternative  is  to  have  an  initially  stoichiometric 
mixture,  in  which  case  the  mixture  would  be  fuel-rich  at  the  end  of 
the  injection  process.   This  did  not  significantly  reduce  the  nitric 
oxide  formation  unless  large  quantities  (more  than  25%  of  the  primary 
fuel)  of  the  mixture  were  injected,  in  \/hich  case  the  carbon  monoxide 
emissions  rose  sharply.   Tliis  is  also  bound  to  have  an  adverse  effect 
on  performance  by  affecting  the  flammability  of  the  fuel-air  mixture. 

It  is  thus  seen  that  this  type  of  an  injection  process  is 
not  particularly  effective  in  reducing  the  emissions,  the  reduction 
of  one  of  the  pollutants  (NO)  adversely  affecting  another  pollutant 
(CO).   This  is  in  general  true  of  most  of  the  techniques  to  control 
the  emission  of  the  oxides  of  nitrogen.   The  difficulty  can  be  over- 
come by  reducing  the  oxides  of  nitrogen  to  the  specified  levels 
by  these  control  techniques  and  later  on  providing  catalytic  treat- 
ment of  the  exhaust  gases.   The  latter  is  now  considerably  simplified 
and  involves  the  conversion  of  carbon  monoxide  to  carbon  dioxide  by 
an  amount  v/hich  would  reduce  the  former  to  the  specified  levels. 
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CHAPTER  VII 
CONCLUSIONS  AND  RECOMt-IENDATIONS 

A  theoretical  analysis  has  been  presented  for  the  study  of 
the  thermodynamic  processes  occurring  in  an  Otto  cycle  engine.   A 
detailed  investigation  of  the  combustion  process  has  been  carried 
out  by  analyzing  the  propagation  of  a  flarae  front  into  a  mixture 
of  homogeneous  combustible  gases.   The  solution  technique  presented 
has  been  demonstrated  for  the  case  of  a  methane-air  flarae  and  a 
description  of  the  structure  of  the  flame  zone  has  been  obtained  in 
considerable  detail.   The  technique  may  be  used  to  predict  the 
structure  of  other  flames  depending  on  the  availability  of  chemical 
kinetics  data.   Alternatively,  the  solution  method  provides  a 
powerful  tool  for  obtaining  information  on  the  reaction  mechanism 
and  chemical  kinetics  data  from  experimentally  determined  flame 
structures . 

The  processes  leading  to  the  formation  and  later  partial 
decomposition  of  two  of  the  major  pollutants  from  a  spark-ignition 
engine,  namely  carbon  monoxide  and  nitric  oxide,  have  also  been 
analyzed  in  detail.   The  formation  of  the  pollutants  has  been 
considered  on  the  basis  of  a  finite  rate  combustion  process.   The 
decomposition  process  during  the  expansion  stroke  has  been  studied 
on  the  basis  of  a  non-equilibrium  analysis.   It  has  been  demonstrated 
that  the  assumption  of  chemical  equilibrium  throughout  the  cycle  does 
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not  portray  the  true  state  of  affairs  with  respect  to  the  prediction 
of  the  emission  levels. 

The  problem  under  study  involves  a  consideration  of  several 
areas,  e.g.,  fluid  mechanics,  thermodynamics,  chemical  kinetics, 
transport  processes,  computational  techniques,  etc.   In  this  work, 
a  lack  of  reliable  chemical  kinetics  data  was  felt  acutely  and  much 
work  remains  to  be  done  in  this  area  so  as  to  obtain  more  realistic 
theoretical  predictions.   The  data  so  far  reported  are  for  relatively 
low  tem.perature  ranges  and  future  work  should  concern  temperatures 
encountered  in  actual  internal  combustion  engines  or  temperatures 
approaching  such  values. 

The  same  holds  true  for  the  prediction  of  transport  properties, 
especially  of  free  radicals.   Tne  lack  of  experimental  data  makes  it 
difficult  or  even  impossible  to  evaluate  the  accuracy  of  the 
theoretical  predictions  which  again  are  quite  scanty.   An  accurate 
determination  of  the  transport  properties  would  enable  the  determina- 
tion of  reaction  mechanisms  and  reaction  rate  constants  from 
experimentally  m.easured  flame  structure  data.   A  typical  procedure 
would  be  to  start  with  a  postulated  mechanism,  obtain  the  flame 
structure  in  detail  by  theoretical  calculations,  and  compare  this 
with  the  experinientally  measured  values;  the  process  being  repeated 
until  satisfactory  agreement  betv7een  the  tvra  is  obtained.   VJith  this 
in  view,  a  fairly  general  computational  scheme  has  been  presented 
for  the  theoretical  prediction  of  the  flame  structure. 

A  consideration  of  the  processes  occurring  in  the  exhaust 
pipe  would  lead  to  theoretical  predictions  which  can  be  more  closely 
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matched  to  the  experimentally  measured  values.   The  importance  of 
the  chemical  reactions  in  effecting  further  composition  changes  of 
the  gases  needs  further  study.   The  study  would  be  useful  in  the 
evaluation  of  post-exhaust  treatment  of  gases  as  with  secondary 
air  injection  into  the  exhaust  manifolds  and  the  use  of  catalytic 
converters , 

The  phenomenon  of  heat  transfer  during  the  combustion  and 
the  expansion  processes  requires  further  study.   Due  to  the  complex 
nature  of  the  processes  this  would  represent  a  separate  problem  by 
itself,  but  one  worth  looking  into.   The  analysis  of  this  problem 
would  also  shed  some  light  on  the  structure  of  the  quench  zone,  an 
important  source  of  the  unburnt  hydrocarbons  in  the  exhaust. 

Finally,  the  fluid  mechanical  aspects  of  the  gas  flow  into 
and  out  of  the  engine  cylinders  deserve  some  mention.   Some  important 
factors  which  require  further  study  are:   (a)  distribution  of  the 
fuel-air  mixture  in  the  various  cylinders  leading  to  non-uniform 
air-fuel  ratios  and  intake  manifold  design;  (b)  the  fluid  mechanics 
of  the  exhaust  process  to  determine  ivfhat  percentage  of  the  unburnt 
and  partially  burnt  hydrocarbons  from  the  quench  zone  are  actually 
expelled  from  the  engine  cylinders;  the  effect  of  valve  overlap; 
and  (c)  the  study  of  the  mixing  between  the  quench  zone  gases  and 
the  burnt  gases  after  flam.e  propagation  leading  to  partial  oxidation 
of  the  hydrocarbons  in  the  initial  quench  zone,  i.e.,  the  effect  of 
the  level  of  turbulence  on  the  mixing  and  reactions  in  the  post-flame 
gases. 


EQUILIBRIUII   CAIXULATIONS 


Introduction 


Several  papers  dealing  vjith  the  cop.putation  of  chemical 
equilibria  have  appeared  in  the  literature.   However,  most  of  the 
methods  suggested  are  simply  minor  modifications  of  the  previous 
methods,  the  improvements  being  more  related  to  the  numerical 
analysis  than  to  the  principle  of  the  m.ethod  of  solution. 

Tnere  are  two  main  methods  available  for  the  com.putation  of 
chemical  equilibria.   The  first  one  involves  the  minimization  of 
Gibbs'  free  energy.   The  second  method  involves  the  solution  of 
the  set  of  non-linear  algebraic  equations  obtained  from  mass 
balance  and  the  law  of  mass  action.   Several  numerical  techniques 
for  the  solution  of  such  equations  have  been  reported.   The  method 
used  here  is  the  latter  one,  and  is  basically  the  one  due  to  Kandiner 
and  Brinkley  [1]  with  some  modifications  [2]. 

The  principle  of  this  method  is  to  (a)  select  a  certain 
nuiii)er  of  species  (referred  to  as  'independent  components'); 
(b)  estimate  their  equilibrium  concentrations;  and  (c)  express  the 
concentrations  of  the  rest  of  the  species  (referred  to  as  'derived' 
constituents  or  species)  in  terms  of  the  concentrations  of  the 
independent  components  through  the  use  of  chemical  equilibrium 
constants.   Tlie  concentrations  of  the  derived  species  are  then  used 
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as  corrections  in  the  atom  balance  equations  to  obtain  better 
estimates  of  the  concentrations  of  the  independent  components;  the 
process  is  repeated  to  the  desired  degree  of  convergence. 

The  problem  here  is  thus  to  find  the  equilibrium  concentra- 
tions of  the  various  chemical  species  (assumed  16  in  number:   C0„ , 
CO,  N^,  H^O,  H,  OH,  0,  NO,  NO^ ,  N^O,  N,  0^ ,  H2 ,  CH^,  CHO,  CH^)  at  a 
given  temperature  and  pressure  for  given  initial  concentrations  of 
the  reactants. 

For  a  system  at  chemical  equilibrium  comprising  of  N  species 

made  up  of  N  different  types  of  chemical  elements,  the  composition 

of  a  chemical  species  S.  may  be  expressed  in  the  form: 

N 
c 

S.  =  y  a.   A      (i=l,2, N)      (1) 

1    '"^  ic   c 
c=l 

Since  there  are  N  elem.ents  present  in  the  svstem,  there  are 
c 

at  the  m.ost  N  linearly  independent  vectors  a.  .   Tlie  (N-N  )  linearly 
c  ^  1  c 

dependent  vectors  are  linear  combinations  of  these: 

N 
c 

ya.v.=a..  (2) 

^,   IC  IC      11 

c=l    -^ 

corresponding  to  the  chemical  relations: 

K 
c 

y  V.  S  =  S.  (3) 

'^^  ic  c     1 
c=l 

Kandiner  and  Brinkley  [1]  use  the  term  independent  components  to 
describe  the  species  corresponding  to  the  linearly  independent 
vectors  a. . 
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ITie  first  step  in  this  method  involves  the  determination  of 
the  number  of  independent  components.   This  number  is  always  less 
than  or  equal  to  the  number  of  chemical  elements  present  in  the 
system.   The  two  are  generally  equal,  an  exception  being  the  case 
when  certain  atoms  appear  only  in  groups  which  go  through  the  chemical 
reactions  intact. 

To  facilitate  numerical  computations,  species  which  are 
expected  to  be  present  in  large  concentrations  at  equilibrium  are 
chosen  as  the  independent  components.   This  choice  is  subject  to  the 
following  conditions: 

(i)   The  species  so  chosen  must  be  stoichiometrically  in- 
dependent of  one  another,  e.g.,  the  species  CO  ,  CO,  0   cannot  all 
be  selected  as  the  independent  components  since  CO  =  CO  +  ~  0^;    any 
two  of  the  three  may  be  selected. 

(ii)   Tne  group  of  species  selected  as  independent  should 
contain  all  the  chemical  elements  present  in  the  system. 

(iii)   If  a  condensed  phase  is  present,  it  must  be  chosen  as 
an  independent  component  (e.g.,  solid  carbon). 

In  this  case,  the  independent  com.ponents  are  four  in  number 
(equal  to  the  number  of  chemical  elements  present — C,  0,  H,  N) ;  these 
are  chosen  to  be  CO   CO,  N^ ,  HO  since  these  are  expected  to  be 
present  in  the  largest  concentrations  in  the  final  equilibrium 
mixture.   Also  all  of  the  above  conditions  are  satisfied.   The 
remaining  species  in  the  system  are  referred  to  as  the  derived 
constituents  or  derived  species.   It  is  assumed  that  solid  carbon 
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is  absent  from  the  equilibrium  mixture.   The  presence  of  condensed 
phases  is  discussed  below. 

The  next  step  involves  the  writing  of  the  equilibrium 
relations  for  each  one  of  the  derived  species,  such  that  in  each 
relation  only  the  corresponding  derived  species  and  any  (or  all) 
of  the  independent  components  are  present.   These  relations  are: 

H:  i  "2°  "^  i  ^°  "  i  ^^2  "^  ^ 

OH:  i  "2°  "^  I  "^^2  =  f-  CO  +  OH 

0:  CO  =  CO  +  0 

NO:  ~   N  +  CO  =  CO  +  NO 

NO^:        i  ^2  "*"  ^^°2  "  ^*^°  "^  ^'°2 
N^O:        N^  +  CO2  =  N^O  +  CO 

N:  i  N2  =  N 

0  :  200^  =  2C0  +  0^ 

CH^:  I.5H2O  +  3, SCO  =  CH^  +  2.5CO^ 

CHO:  0.5H  0  +  l.SCO  =  CHO  +  0.5CO 

CH,  :        2H^0  +  4C0  =  CH,  +  3C0,, 
4  2  4      2 

Additional  relations  to  be  used  when  C  H  „,  CH  OH  and  C  H^OH 
are  considered: 

CgH^g:      9H2O  +  25CO  =  CgH^^  +  1700^ 

CH  OH:       3C0  +  2H  0  =  200^  +  CH  OH 

C  H  OH:     6C0  +  3H  0  =  4C0  +  C  H  OH 
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For  each  one  of  these  equilibrium  relations,  the  law  of  mass 
action  is  of  the  form: 


PJ 


I  V,   -  1 


p  ^c   -^ 


n  X 


JC 


or, 


X.  =  K  .  p 
3  PJ 


I  V.   -  1 

C   '' 


V. 

n  X  J'^ 


(4) 


where 


log  „  K  .  =  )]  V.   log^^  K  -  -  log,„  K  ^. 
10   PJ      jc   ^10   pf     ^10   pfj 


Values  of  log.Q  K   are  obtained  from  the  JANAF  Tables. 
Subscript  j  above  refers  to  the  derived  constituents; 


J  =  (\  +  1)  ,   (\  +  2)  , 


Subscript  c  refers  to  the  independent  components; 


c  =  1  ^ 


This   notation   is    used   throughout    the    remainder   of    this    analysis 
'Ihe    atom  balance   equations    are: 


Ta.   n.+Ta     n     =B 
.      nei         ''cec  e 

J  c 


(5) 


for  all  elements  e,  e  =  1,2, N  ;  and 


B  =  total  number  of  atom.s  of  element  e 
e 


which  is  known  from  the  initial  system  composition. 
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An  arbitrary  reference  element  e  =  r  is  chosen  and  for  this 
element , 


ya.n.+yan=B  (6) 

.   ir  1   ^   cr  c    r 
-J   J   J    c 


From  equations  (5)  and  (6), 


Tn  (a.      -Ra.)+yn   (a   -Ra  )=0      (7) 
.      1        TB  e  ir    ^   c   ce    e  cr 


where   R  =  B  /B 
e    e   r 

Dividing  by  n  =  )]  n   +  )]  n^  , 


y  X.  (a.   -  R  a.  )  +  y  X   (a   -  R  a   )  =  0      (8) 
^   1    ie    e  ir    ^   c    ce    e  cr 
3  c 


Also  the  sum  of  the  mole  fractions  of  the  independent  and  derived 
components  equals  unity: 


I  X.  +  I  X  =  1  (9) 

.1    ^   c 
J      c 


The  problem,  is  now  reduced  to  the  simultaneous  solution  of 
equations  (A),  (8)  and  (9),  which  together  comprise  a  set  of  N 
simultaneous  algebraic  equations  in  the  unknowns  x.  (i=l,2, N) , 

The  Newton-Raphson  linearization  method  is  used  to  obtain  the 
solution  of  these  equations.   The  method  is  outlined  here;  details 
may  be  found  in  Ref.  2.   Tliis  method  involves  estimating  the  values 
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of  the  mole  fractions  of  the  independent  components;  the  estimated 
values  for  the  derived  constituents  are  then  obtained  using  equations 
(4).   Denoting  the  estimated  solution  by  (y.,y  ),  equations  (8)  and  (9) 
are  in  general  not  satisfied  by  (y.,y  ),  and  so, 


I   y.  (a.,  -  R.a.    )   +  I  y      (a  ,  -  R,  a   )  ^ 
.2        2^  K  jr    ^  •'c   ck    k  or 


(10) 


=  Fj_   (say) , 


k  =  1,2, (N^-  1) 


and     I   Y.   +  )   y     -It^O 
J      c 

=  7^      (say)  ,   k  =  N^ 


(11) 


These  N  equations  are  written  as  one  general  equation  of  the  form: 


^1  =Ia-iy-+Ia,y 
k    .   ik  1    ^   ck^c 
J         c 


(12) 


wnere 


a.,  =  a.,  -  R,  a. 
jk    jk    k  jr 


a  ,  =0 
ck    cl' 


\^ 


cr 


[   k=  1,2,, 


(\-  1) 


(13) 


a.,  =  1 


-ck=l 


\  =  1  +  \ 


k=  N 


(14) 


Also  from  equations  (4), 


y.  =  K.  n  y 
3     J  c   c 


V. 


(15) 
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V.   -  1 

where         K.  =  K  .  p 
J     PJ 


Considering  the  y   to  be  independent  variables,  the  functions 
F  may  be  expanded  in  a  Taylor  series,  and  finally  the  following 


relations  may  be  obtained  [2]: 


0  =  Fv  +  ^ 


—           V.  — 

a  ,  +  y  a.,  -^^  •  y. 

ck    .   jk  y  ^1 
J      c 


6y  (1&) 

c 


This  is  a  set  of  N  equations  for  the  correction  terms  Sy    . 

c  c 

Tliese  equations  are  solved  using  Gauss-Jordan  elimination  with 
maximum  pivoting.   The  new  estimates  for  the  y   are  then: 


(y  )    =  y  +  6y  (17) 

c  new   ^  c     c 


The  process  is  repeated  till  the  mole  fractions  of  all  the 
species  converge  to  within  a  predetermined  accuracy. 


Convergence 


Tiie  amplitude  of  each  successive  error  is,  roughly  speaking, 
proportional  to  the  square  of  the  amplitude  of  the  preceding  error. 
The  main  difficulty  v;ith  this  method,  apart  from  the  need  to  cal- 
culate and  invert  a  matrix  at  each  step,  is  that  the  starting  estimate 
should  be  a  rather  good  one,  if  the  method  is  to  converge.   A  poor 
starting  estimate  largely  invalidates  the  analysis  since  third  and 
higher  order  terms  are  no  longer  negligible.   If  the  initial 
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estimates  are  poor,  such  that  some  of  the  mole  fractions  have  been 
largely  overestimated,  the  corrections  may  be  so  largely  negative, 
that  negative  mole  fractions  are  calculated.   The  introduction  of 
logarithmic  variables  overcomes  this  problem.   Since  this  method 
involves  a  first  order  approximation,  one  can  use  as  an  approximation 
for  6y  /y   the  function  6£ny  which  is  equal  to  6y  /y   to  first  order. 
Then, 


(£ny  )    =  £ny  +  5£ny  (18) 

c  na\<i  c       c 


(v  )    =  y  exp  (6£ny  )   i  0  (19) 

c  new   -^  c  c   ^ 


The  Newton-Raphson  method  converges  rapidly,  if  it  does 
converge  at  all.   However,  depending  on  the  initial  estimate,  the 
iteration  may  diverge  or  there  may  not  be  even  a  unique  solution. 


Vfnen  the  determinant  of  the  function  F 


'^ck 


V. 

a  ,  +  y  a.,  -^^  y . 
ck   V   |k  y   ■'2 


given  by 
vanishes  or  is  close  to  zero,  the  matrix 


becomes  singular  or  near  singular  and  the  solution  would  exhibit 
divergence , 

Kandiner  and  Brinkley  [1]  have  used  m.ole  numbers  in  their 
calculations  instead  of  mole  fractions  used  here.  There  are  some 
difficulties  involved  in  the  use  of  mole  numbers: 

(i)   An  extra  equation  I  n .  +  )^  n   =  n  has  to  be  solved  at 
each  stage  of  the  com.putation . 

(11)   Tne  initial  estimates  have  to  obey  the  atom  balance 
equations;  this  needs  the  inversion  of  a  m.atrix. 
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(iii)   Negative  mole  numbers  are  sometimes  calculated  because 
of  overcorrection  caused  by  poor  estimates. 

(iv)   Mole  fractions  are  easier  to  estimate  than  mole  numbers. 

Inclusion  of  condensed  phases 

The  original  treatment  of  Kandiner  and  Brinkley  [1]  was  a 
very  general  one  in  that  it  included  consideration  of  heterogeneous 
systems.   In  the  case  of  the  assumed  presence  of  condensed  phases, 
the  procedure  is  basically  identical,  in  that  the  equations  are  all 
applicable  in  each  phase.   Each  solid  species  has  to  be  selected 
as  one  of  the  independent  components.   The  maximum  number  of  condensed 
phases  present,  in  the  case  where  a  gas  phase  is  present,  is  equal  to 
(N  -  1) ,  in  which  case  all  solids  and  one  gaseous  constituent  are  to 
be  selected  as  the  independent  components ,  and  the  problem  is  very 
much  simplified.   Tnis  restricts  the  number  of  condensed  species  that 
may  occur,  but  does  not  restrict  the  total  number  of  condensed  species 
that  can  be  considered  for  the  solution.   Tlie  main  problem  in  dealing 
with  condensed  phases  is  that  an  initial  assumption  has  to  be  made 
as  to  which  phases  are  present.   Tne  correctness  of  this  assumption 
is  determined  only  after  an  iterative  computation  has  converged  by 
checking  whether  the  vapor  pressure  of  the  particular  species  in  question 
is  equal  to  or  less  than  the  saturation  vapor  pressure  at  that  temper- 
ature.  If  wrong,  another  assumption  must  be  made  and  the  computation 
must  be  repeated. 
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NOMENCLATURE 


a.    =  vector  of  a. .    for  given  1 

=  ^^l'^2' X^ 

A    =  chemical  element  c 
c 

B    =  number  of  atoms  of  chemical  element  e  (e=l,2, N  ) 

K.    =  equilibrium  constant  in  terms  of  mole  fractions 


K  . 
PJ 


equilibrium  constant  in  terms  of  partial  pressures 
N    =  number  of  chemical  species 


N    =  number  of  chemical  elements 


n 


mole   num.ber   of   derived  species   j    (j=N     +  1, N) 


J  r  u   v^    ^ 


n 
c 


mole  number  of  independent  component  c  (c=l,2, N  ) 


p  =   pressure  in  atmospheres 

R  =   B  /B   (e=l,2, N  ) 

e       e   r     '  '       c 

S.  =   chemical  species  i 

V 


ij 


stoichiometric  coefficient  of  independent  component  j  in 
the  equation  for  the  formation  of  derived  species  i 


X.  =  mole  fraction  of  derived  species  j  (j=N  +  1, N) 

-J  ^ 

X  =   mole  fraction  of  indeoendent  component  c  (c=l,2 1<  ) 

c  '  '^  '  '     c 

y.  =   estimated  value  for  y. .     (j=l,2, N) 

Ct.  -■  number  of  atom.s  of  element  c  in  soecies  i 
ic 


Subscripts: 

C    =   independent  component 

j     =   derived  species 

r    =   reference  element 
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APPENDIX  II 
THE  DETERI-IINATION  OF  TRANSPORT  PROPERTIES 

The  conservation  equations  for  the  analysis  of  the  flame  zone 
contain  terms  associated  with  the  fluxes  of  mass,  energy  and  momentum 
due  to  gradients  in  the  m.acroscopic  properties.   The  coefficients 
of  diffusion,  thermal  conductivity  and  viscosity  which  appear  in 
the  conservation  equations  m.ust  first  be  determined  in  order  to  solve 
these  equations.   All  these  three  phenomena  are  physically  similar 
in  the  sense  that  they  involve  the  transport  of  some  physical  property 
(mass,  heat,  momentum)  through  the  medium  under  consideration. 

In  this  appendix  a  review  of  the  status  of  the  transport 
properties,  as  applicable  to  the  problem  under  study,  is  made. 
Details  of  the  derivations  of  the  various  equations  are  not 
considered  and  may  be  found  in  the  appropriate  literature.   The 
main  emphasis  is  on  the  results,  their  limitations,  the  various 
assumptions  involved  and  how  far  these  assumptions  are  met  in  this 
problem.   Dixon-Lewis  [1]  has  presented  an  excellent  reviev?  of  the 
literature  on  transport  phenomena  as  related  to  flame  studies. 

Detailed  experimental  data  on  the  transport  coefficients 
of  gases  (especially  the  diffusion  coefficient)  even  at  room 
temperature  are  not  available  and  data  at  temperatures  encountered 
in  flames,  or  even  at  temperatures  approaching  such  values,  are 
virtually  non-existent.   Very  few  measurements  at  high  temperatures 
have  been  made  and  even  these  are  not  very  accurate.   The  reason  is 
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that  measurements  even  at  ordinary  temperatures  are  difficult  and 
high  temperatures  make  the  problem  worse .   Experimental  data  on 
the  transport  properties  in  multi-component  systems,  as  encountered 
in  flame  studies,  are  almost  completely  lacking,  with  only  a  few 
approximate  measurements  having  been  reported. 

This  lack  or  scarcity  of  experimental  data  has  made 
theoretical  prediction  of  transport  coefficients  (especially  at 
high  temperatures)  the  best  available  estimates.   Among  the  different 
theoretical  methods  one  of  the  best  estimates  of  the  transport  coef- 
ficients is  that  based  on  the  Chapman-Enskog  theory. 

The  determination  of  the  transport  coefficients  from  the 
rigorous  kinetic  theory  of  gases  is  based  on  the  knoxvf ledge  of  the 
distribution  function  f..  (r,V.  ,  t)  ,*  which  is  obtained  by  solving  the 
Boltzmann  integro-dif ferential  equation.   The  solution  reduces  to 
the  Maxx'/ellian  distribution  for  a  system  at  equilibrium  (i.e.  ,  in  the 
absence  of  concentration,  velocity  and  temperature  gradients). 


i.\    ^   =  n.  (m./27TkT)-''^  exp  (-m.v:/2kT) 


It  is  assumed  in  most  combustion  problems  that  the  gradients 
of  the  various  properties  are  not  so  great  as  to  make  it  impossible 
to  define  a  'simple'  transport  coefficient — i.e. ,  the  conditions  in 
the  problem  are  assumed  to  be  not  too  far  rem.oved  from  equilibrium 
so  that  the  flux  vectors  are  linear  in  their  derivatives.   This  is 


Tnis  function  represents  the  number  of  molecules  of  the  ith 

species  which  at  time  t  lie  in  a  unit  volume  about  the  point  r  and 

which  have  velocities  within  a  unit  range  of  V. . 
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true  when  changes  in  macroscopic  gas  properties  over  the  distance 
of  a  mean  free  path  are  small.   This  is  generally  a  good  assumption 
for  laminar  flame  propagation  but  breaks  down  for  the  case  of  shock 
waves,  detonations,  etc.,  where  continuum  analysis  ceases  to  be  valid. 

For  near  equilibrium  systems  Chapman  and  Enskog  have  used  a 
series  expansion  of  the  distribution  function  about  the  Max^^7ellian 
distribution,  and  solved  the  Boltzmann  equation  by  a  perturbation 
method.   This  finally  leads  to  expressions  for  the  transport 
coefficients  in  terms  of  a  set  of  integrals  ^  '^\    called  the 
collision  integrals  (1  =  number  of  molecules  which  collide 
simultaneously;  s  =  number  of  molecules  affecting  the  behavior  of 
the  collision;  s  i  1).   Tnese  integrals  depend  on  the  dynamics  of 
the  intermolecular  collisions  and  hence  on  the  intermolecular  forces. 
The  nature  of  these  forces  is  not  very  well  known,  especially  for 
complex  (and/or  polar)  molecules  and  this  limits  the  applicability 
of  the  theoretical  predictions. 

For  non-polar  miolecules  the  empirical  Lennard- Jones  (6-12) 
intermolecular  potential  energy  function  is  widely  used.   This  is: 


(J)(r)  =  4c 


12 

a 

r 


(1) 


The  parameters  O   and  e  have  the  dimensions  of  length  and 
energy  respectively.   These  are  constants  for  particular  chemical 
species.   The  physical  significance  of  O   and  £  is: 

(1)   r  =  a  denotes  the  closest  approach  of  two  molecules  which  collide 
witl'i  zero  initial  relative  kinetic  energy  since  <^(r)   =  0  at  r  =  a- 
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hence  the  name  of  'collision  diameter'  for  O  is  used  extensively. 
(2)  e  denotes  the  maximum  energy  of  attraction  between  a  pair  of 
molecules  and  this  occurs  at  a  separation  of  r  =  2"*"   'a. 

To  determine  the  viscosity  and  thermal  conductivity  it  is 
necessary  to  have  values  of  O   and  £  and  for  the  diffusion  coefficient 
values  of  O^^    and  e^^ ,  where  o^      and  e^     are  the  Lennard-Jones  parameters 
for  a  pair  of  unlike  molecules  i  and  j . 

The  Lennard-Jones  potential  is  quite  often  a  good  approximation 
even  for  molecules  which  are  not  spherical  or  non-polar.   Since  the 
coefficient  of  viscosity  is  not  strongly  affected  by  the  internal 
degrees  of  freedom  it  is  possible  to  find  the  Lennard-Jones  potential 
parameters  by  fits  to  experimental  data  for  these  molecules.   Tnus 
this  potential  function  is  used  for  non-spherical  and  polar  molecules 
where  the  experimental  data  warrant  its  use.   However,  for  certain  gases, 
the  experim.ental  data  cannot  be  fitted  to  this  potential  and  recourse 
must  be  taken  to  other  potential  functions.   Notable  among  such 
molecules  and  of  interest  in  combustion  studies  are:   water  which  is 
strongly  polar,  and  long  molecules  such  as  higher  hydrocarbons.   For 

example,  x/ith  water-air  mixtures  (considered  as  binary)  the  experimental 

1 .  89 
data  follow  a  T  '    dependence  (T  =  temperature),  while  the  Lennard-Jones 

potential  function  predicts  a  T  '    dependence,  causing  serious  errors 
at  high  temperatures. 

The  Buckingham  potential  (6-exp)  has  also  been  used  to  a  certain 
extent  to  describe  intermolecular  forces.  It  has  three  adjustable  para- 
meters instead  of  the  two  in  the  Lennard-Jones  potential  and  should 
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therefore  in  principle  allow  better  fit  of  experimental  data.   The 
advantages  however  are  almost  negligible  [2,3].   Therefore,  this 
potential  function  xjill  not  be  discussed  any  further. 

Diffusion  coefficient.   Firstly,  brief  mention  v/ill  be  made 
of  the  experimental  data  which  are  available.   The  first  high  temper- 
ature measurement  of  diffusion  coefficient  was  reported  by  Klibanov 
et  al.  [4].   Their  measurements  for  the  diffusion  coefficient  of 
CO^-air  and  H^O-air  mixtures  (considered  as  binary  mixtures)  X'7ere 
made  at  temperatures  up  to  about  1500 °K.   The  data  are,  however, 
highly  imprecise,  as  to  be  of  no  practical  significance  from  a 
quantitative  standpoint.   Other  high  temperature  measurements  are 
those  reported  by  Walker  and  Westenberg  [5-8]  and  Westenberg  and 
Frazier  [9]  for  gases  of  interest  in  the  study  of  methane-oxygen 
flam.es,  over  the  temperature  range  300-1150°K,  except  for  some  gas 
pairs  V7hich  underwent  chemical  reactions  at  lower  temperatures. 
Ember,  Ferron  and  Wohl  [10]  have  reported  self-diffusion  data  for 
CO2  at  temperatures  up  to  1700°K.   Other  published  data  [11,12] 
have  been  on  gases  not  of  interest  in  combustion  problems  and  even 
these  are  not  at  very  high  temperatures,  approximately  eOCK.   The 
precision  of  experimental  data  on  the  diffusion  coefficient  varies; 
the  recent  data  are  often  accurate  to  within  1%,  while  the  accuracy  of 
the  older  data  is  generally  within  10%.   The  ex-perimental  errors 
generally  give  values  of  the  diffusion  coefficient  which  are  higher 
than  the  correct  value. 

The  theoretical  prediction  of  the  diffusion  coefficient  based 
on  the  Chapman-Enskog  theory  and  the  Lennard- Jones  potential  requires 
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a  knowledge  of  the  parameters  O..    and  e . . .   In  principle,  the  best 
way  of  obtaining  these  parameters  is  by  an  empirical  fit  to  accurate 
measurements  of  the  diffusion  coefficient  over  a  wide  range  of 
temperatures.   Due  to  the  absence  of  suitable  experimental  data 
this  has  been  done  in  only  a  very  fe\<r   cases  (for  example,  see 
Refs,  3,  13).   For  non-polar,  non-reacting  molecule  pairs  fairly 
satisfactory  estimates  can  be  obtained  [14]  by  the  usual  empirical 
combining  laws : 


a,.  =i(a.  +a.)  (2) 


where  O^,   o        e^,    e      are  obtained  by  empirical  fit  to  the  much  m.ore 
abundant  data  on  viscosity  of  the  individual  species  or  from  second 
virial  coefficient  data.   However,  it  should  be  noted  that  the 
empirical  fits  to  one  property  of  a  pure  gas  need  not  necessarily 
be  the  best  one  for  determining  another  property  of  a  mixture  of 
that  gas  with  other  gases.   For  these  calculations,  it  is  generally 
better  to  use  viscosity  data,  if  available,  and  take  the  values  of 
a  and  £  from  the  same  source  (if  possible)  rather  than  use  second 
virial  coefficient  data.   The  restriction  that  the  molecule  pair  be 
non-reacting  has  been  mentioned  above.   Even  for  reacting  gases  the 
number  of  collisions  resulting  in  reaction  is  very  sm^all  compared  to 
the  total  number  of  collisions,  so  these  relations  may  be  used. 
However,  when  the  reactions  are  very  rapid  (half-life  of  approxim.ately 
10  microseconds)  complications  are  introduced  since  a  much  larger 
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fraction  of  the  total  number  of  collisions  results  in  chemical  reaction. 
These  complications  are  not  considered  here  and  the  simple  combining 
rules  are  assumed  to  hold  even  though  rapid  chemical  reactions  are 
occurring  in  the  flame  zone. 

The  viscosity  data  generally  give  lower  values  of  the 
diffusion  coefficient  as  compared  to  the  values  obtained  by  best 
fits  to  experimental  diffusion  data  and  extrapolated  to  the  required 
conditions.   At  high  temperatures  (1000  degrees  K,  or  higher),  this 
error  may  amount  to  as  much  as  20-30%  for  some  cases  (but  in  many 
cases  they  are  more  precise),  while  at  lower  temperatures  the 
accuracy  is  much  better,  of  the  order  of  5%,  which  is  generally 
within  the  limits  of  experimental  error.   Due  to  lack  of  experimental 
data,  the  high  temperature  predictions  have  not  been  adequately  tested. 

I'Jhen  determining  the  parameters  O   and  e   by  empirical  fits,  it 
should  be  noted  that  the  collision  integrals  are  little  affected  by 
the  uncertainties  in  the  values  of  e.      However,  the  parameter  O 
occurs  as  a  squared  quantity  in  the  evaluation  of  the  transport 
properties.   Low  values  of  a  can  be  offset  by  high  values  of  E  and 
vi  ce  ve  rs  a . 

For  nonpolar-polar  gas  pairs  the  Lennard-Jones  potential 
is  still  used,  since  apparently  the  polarity  of  one  of  the  species 
does  not  strongly  influence  the  molecular  interaction  and  approximate 
spherical  symm.etry  nay  still  be  assumed.   Hov/ever,  the  above  combining 
laws  have  to  be  modified: 
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a   =  ^  (a  +  a  )  f  ^/^  (4) 

np   2   n    p    c  ^    ' 


e       =  (e  £  )^/^  f^  (5) 

np     n  p      c  ^  ^ 


f  =  1  +  (l/v^)(a  t*/a-^)/£  /£ 
c  ^  n  p'  n'   p'  n 


where   a  =  the  polarizability  of  the  non-polar  molecule, 

y  =  dipole  moment  of  the  polar  gas, 

n,p  =  subscripts  for  non-polar  and  polar  molecules, 
respectively , 


(6) 


and 


t*  =  1//8  •  (y^/£  a^) 
P  P  P  P 


(7) 


For  polar-polar  gas  pairs  the  Stockmayer  potential  is  used 
for  the  estimation  of  the  diffusion  coefficient.   This  potential 
energy  function  is  the  sum  of  the  Lennard-Jones  function  and  an 
angle-dependent  function  proportional  to  the  inverse  third  po^v'er  of 
intermolecular  separation.   This  function  takes  care  of  the  electro- 
static interaction  between  two  dipoles.   The  theory  of  transport 
phenomena  in  polar  gases  has  been  developed  by  Mason  and  Ilonchick 
[15,16].   The  final  results  are  g'iv^  in  terms  of  the  collision 

integral  fi  which  is  tabulated  [16]  as  a  function  of  kT/£   and 
u  P 

(6)    where: 
max 


(6)    =  i  y^/a^  e  =  /It*  (8) 

max   ^   P   P   P       P 


ITie  combining  rules  for  £.  .  and  O .  .    are  the  same  as  those 

given  by  equations  (2)  and  (3)  and  for  (6..)    : 

ij  max 
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(5..)    =  ^  M.y./e..a^  =  [(6.)    (6.)   ]^''^[{o .o  )^'^ /o .    ]'^ 
xj  max   2  I'^j   ij  13      1  max   j  ^max    ^^  i  j^     ij 

(9) 

Due  to  the  virtual  non-existence  of  data  for  polar  gas  pairs 
it  is  difficult  to  predict  the  accuracy  of  this  theory.   The  potential 
functions  described  above  are  reasonably  adequate  for  simple  molecules. 
However,  they  cannot  be  used  to  describe  interactions  involving  long 
molecules,  free  radicals,  ions,  and  molecules  in  excited  states. 

Thermal  conductivity.   The  amount  of  experimental  data  available 
for  thermal  conductivity  is  considerably  more  than  that  for  diffusion. 
An  extensive  compilation  of  the  data  may  be  found  in  Ref,  17.   The 
data  cover  a  fair  range  of  temperature  and  the  em.pirical  fits  to  the 
data  are  reliable  to  x^^ithin  2-5%  in  the  temperature  range  for  which 
the  data  were  gathered.   The  experimental  determination  of  thermal 
conductivity,  especially  at  high  tem.peratures ,  is  complicated  by 
other  factors  such  as  radiation.   An  exact  consideration  of  radiation 
and  other  losses  is  often  not  possible  and  this  has  led  to  a  wide 
scatter  in  the  experimental  data  gathered  by  different  Investigators. 
Scatter  of  as  much  as  10%  at  temperatures  of  approximately  1000 
degrees  K  with  higher  scatter  at  higher  temperatures  is  not 
uncommon  for  some  of  the  experimental  data  [18]. 

V/ith  viscosity  parameters  O   and  £  the  Chapman-Enskog  theory 
gives  fairly  good  agreement  x^^ith  experimental  values  for  monatomic 
gases  ever  the  m.odest  temperature  ranges  for  which  reliable 
experimental  data  are  available. 

Difficulties  arise  for  polyatomic  molecules,  due  to  the 
presence  of  internal  degrees  of  freedom.   Collisions  between  such 
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molecules  are  inelastic  in  the  sense  that  kinetic  energy  is  no  longer 
conserved,  but  mass  and  momentum  are.   Thus  the  thermal  conductivity 
is  the  only  transport  coefficient  to  be  significantly  affected  by 
the  presence  of  internal  degrees  of  freedom.   For  the  coefficients 
of  diffusion  and  viscosity,  the  theory  of  monatomic  gases  may  be 
applied  even  in  the  case  of  polyatomic  gases  when  the  shape  of  the 
molecules  is  not  too  different  from  spherical,  since  the  assumption 
of  a  spherically  symmetric  potential  of  interaction  has  been  made. 
To  account  for  the  presence  of  the  internal  degrees  of  freedom  the 
semi-empirical  correction  proposed  by  Eucken  is  widely  used: 


^  =  ^  FT  R-  +  4 
^  u 


(10) 


The  modified  Eucken  correction  is  also  used  sometimes.   It 


xs; 


R 


1.32  -^  +  1.77 


(11) 


The  Eucken  correction  is  derived  on  the  assumption  that  the 
distribution  of  species  among  the  internal  energy  states  is  that 
which  is  characteristic  of  local  thermodynamic  equilibrium  and  that 
the  species  in  the  various  quantum  states  have  the  same  diffusion 
coefficient.   In  very  hot  flames  (approximately  AOOO^K,  or  more), 
where  species  may  be  present  in  electronically  excited  states,  both 
these  assumptions  would  be  invalid  and  the  Eucken  correction  would 
not  hold,  and  it  may  even  not  be  possible  to  define  a  unique  thermal 
conductivity  [19]. 
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The  modified  Eucken  correction  yields  values  of  the  thermal 
conductivity  which  are  higher  than  the  corresponding  experimental 
values  while  the  Eucken  correction  yields  values  which  are  lower. 
The  former  is  better  at  higher  temperatures  and  the  latter  at  lower 
temperatures.   Brokaw  [20]  and  Stiel  and  Thodos  [21]  have  given  other 
empirical  correction  formulae  which  represent  a  compromise  between 
the  Eucken  and  the  modified  Eucken  corrections.   Rigorous  analysis 
of  the  phenomenon  requires  knowledge  of  the  translational-rotational 
energy  relaxation  times  and  has  therefore  been  carried  out  only  for 
some  non-polar  gases.   In  the  analysis  carried  out  here,  the  Eucken 
correction  as  given  by  equation  (10)  has  been  used. 

Polarity  causes  more  complications  in  the  determination  of 
thermal  conductivity  than  the  other  transport  properties  since 
energy  transfer  between  polar  molecules  may  involve  internal  (i.e., 
rotational)  modes  and  not  affect  the  translational  energy.   The 
result  is  a  lower  value  of  the  thermal  conductivity.   Mason  and 
Monchick  [22]  have  considered  this  mode  of  energy  transfer  between 
polar  molecules  and  shown  that  the  correction  involved  is  rather 
small  except  for  strongly  polar  molecules  (like  water)  where  the 
difference  betx-zeen  successive  rotational  energy  levels  is  appreciable. 
This  consideration  of  polarity  also  requires  a  knowledge  of  the 
rotational  relaxation  times,  which  is  generally  not  available  for 
most  molecules.   Again,  based  on  these  relaxation  times,  Mason  and 
Monchick  [22]  have  derived  correction  factors  for  polyatomic  gases 
which  do  not  require  the  assumption  of  equilibrium  between  the 
internal  and  translational  energy  modes  (assumption  made  for  the 
Eucken  correction). 
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The  viscosity  of  polar  gases  has  also  been  considered  by 
Mason  and  Monchick  both  for  pure  gases  [16]  and  mixtures  [15]  using 
the  Stockmayer  potential.   They  found  that  the  potential  parameters 
O   and  £  thus  obtained  by  fitting  to  experimental  data  were  not  very 
much  different  from  the  values  obtained  by  fitting  the  data  to  the 
Lennard- Jones  potential  with  the  polarity  neglected. 

Tsederberg  [18]  has  given  a  comparison  of  the  experimental 
data  on  thermal  conductivity  of  di-  and  polyatomic  gases  with  the 
values  calculated  according  to  the  Chapman-Enskog  theory  and  based 
on  the  Lennard-Jones  potential,  and  including  the  Eucken  correction. 
The  comparison  is  made  only  for  temperatures  up  to  300  degrees  C; 
higher  temperatures  were  probably  not  considered  because  of  lack  of 
accurate  experimental  data.   It  is  shown  that  the  agreement  of  the 
ejqjerimental  and  theoretical  values  is  good  for  hydrogen,  oxygen, 
and  carbon  dioxide  with  a  maximum  deviation  of  about  3.7%.   On  the 
other  hand,  for  methane  and  nitric  oxide,  there  is  considerable 
disagreement,  the  maximum  deviation  being  about  10%.   Tne  lack  of 
experimental  data  is  worse  still  for  the  case  of  high  pressures. 
High  pressure  experimental  data  are  available  only  around  room 
tem.peratures  (<  100  °C)  so  that  the  predictions  of  the  high  pressure 
theories  cannot  be  tested. 

Using  the  theoretical  methods,  it  is  possible  to  predict  the 
coefficient  of  viscosity  to  within  2%  for  both  polar  and  non-polar 
gases  as  compared  to  experimental  values.   For  diffusion  coefficients 
the  deviation  is  6-10%  depending  on  the  accuracy  to  v/hich  the 
Lennard-Jones  parameters  a  and  e  have  been  determined,  the  precision 
being  better  if  the  parameters  have  been  determined  from  low 
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temperature  diffusion  data  rather  than  viscosity  data.   Reid  and 
Sherx-7ood  [23,  p.  534]  have  presented  a  comparison  of  the  calculated 
and  experimental  diffusion  coefficients  of  several  gas  pairs.   For 
gas  pairs  of  interest  in  combustion  studies,  some  errors  are  as  high 
as  -17%,  the  average  error  is  about  7%. 

Expressions  for  the  transport  coefficients.   First  approximation 
for  the  viscosity  and  thermal  conductivity  of  a  pure  gas: 


in]      X  10^  =  266.93     "f^  ,  (12) 


[X]   X  10^  =  1989.1     '^   _  (13) 


where   ri  =  viscosity  in  g/cm.sec. 

X   =  thermal  conductivity  in  cal . /cm.sec. °K 
T  =  temperature  in  degrees  K 
T'"  =  reduced  temperature  =  kT/e 
M  =  molecular  weight 
O   =  collision  diam.eter  in  angstroms 
e/k  =  potential  param.eter  in  degrees  K 

For  polyatomic  molecules,  Eucken's  correction  should  be  used 
in  determining  the  thermal  conductivity. 

First  approximation  for  the  diffusion  coefficient  of  a  binary 
mixture  of  species  i  and  j : 
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A-^Qi.    +   M.)/21I.M. 
[D   ]   =  0.0026280  "-  ^|    ^   ^  (U) 

2 

where   D^   =  diffusion  coefficient  in  cm.  /sec. 

p  =  pressure  in  atmospheres 
T  =  temperature  in  degrees  K 
T*.  =  kT/e. . 

^i'^S  ^  molecular  weights  of  species  i  and  j  respectively 

'^,-j>C-;-/k  =  Lennard-Jones  parameters  characteristic  of  i-i 
-■-J   -"-J  -" 

interaction  in  angstroms  and  degrees  K  respectively 
The  higher  approximations  to  the  three  coefficients  are: 


[n]^  =  [n],  f^^'^^  (15) 

n 


The  functions  f  '  ,  f^   ,  and  f^'"^  differ  from  unity  only 

ij 
slightly  and  vary  only  slightly  with  t'\   For  example,  for  the 

(3) 

Lennard-Jones  potential  f '   differs  from  unity  by  less  than  0.8%. 

Tne  first  approximation  to  the  diffusion  coefficient  involves  only 
the  interaction  between  unlike  mLolecules.   Interactions  between  like 
molecules  have  been  neglected  to  a  first  approximation  and  so  there 
is  no  concentration  dependence  term  appearing  in  the  above  expression 
(14).   The  interaction  between  like  molecules  is  considered  in  the 
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second  and  higher  approximations  which  are  complicated  functions  of 
molecular  weights,  mixture  composition  and  collision  integrals  for 
both  like  and  unlike  molecule  interactions.   Experimental  data  to 
verify  this  dependence  on  mixture  composition  are  not  available. 
Explicit  expressions  for  the  higher  order  approximations  (up  to  the 
third)  have  been  obtained  in  terms  of  the  collision  integrals 
[24,  Table  I-P].   These  expressions  show  that  the  dependence  of 
the  diffusion  coefficient  on  the  mixture  composition  is  slight. 
Even  for  gas  pairs  with  very  unlike  masses,  the  higher  order  terms 
amount  to  a  correction  of  a  few  percent.   For  the  Lennard-Jones 
potential,  the  correction  varies  from  1,00  to  1,03  for  most  binary 
gas  pairs  of  interest  in  combustion  studies  and  so  is  almost  always 
neglected. 

Since  the  diffusion  coefficient,  to  a  first  approximation, 
involves  only  the  interaction  betxreen  unlike  molecules,  a  fit  to 
experimental  diffusion  data  is  the  easiest  and  most  accurate  way  of 
determining  the  force  constants  for  interaction  between  pairs  of 
unlike  molecules. 

Approximations  for  the  collision  integrals.   Over  the  range  of 
3  ^  T  ^  200  (generally  encountered  in  flame  work)  ,  the  following 
approximations  may  be  used  [25]: 


P/^'^^  ^   1.12(T''')~°-^^  (18) 


fi^2'2^  ^   1.22(T*)-°-l^  ^^g^ 
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These  are  accurate  to  within  2%  over  the  range  specified. 
Both  the  approximations  break  doxm  rapidly  for  T*  <  3.   Another 
approximation  (slightly  more  accurate)  to  9.^    ">    '   is  [26]: 

1/J^*^^'^^  =  0.697  (1  +  0.323£nT*)  (20) 

Transport  properties  of  mixtures.   The  only  transport  property 
of  mixtures  that  appears  in  the  conservation  equations  is  the  thermal 
conductivity.   For  mass  transport,  only  the  binary  diffusion  coef- 
ficients appear  in  the  diffusion  equations.   As  such  the  discussion 
here  is  limited  to  the  determination  of  the  thermal  conductivity 
of  mixtures. 

The  thermal  conductivity  of  a  gas  mixture  is  not  a  linear 
function  of  composition  in  the  sense  that  it  cannot  be  determined 
by  an  additive  rule  of  the  thermal  conductivities  of  the  individual 
component  gases  multiplied  by  their  respective  mole  or  mass  fractions. 
The  thermal  conductivity  is  proportional  to  the  mean  free  path  of  the 
molecules,  and  in  a  mixture  of  two  or  more  gases  the  mean  free  path 
of  the  molecules  of  each  gas  is  changed  due  to  the  presence  of  the 
molecules  of  all  the  other  gases.   For  mixtures  of  polar  gases,  if 
the  different  types  of  molecules  differ  greatly  in  polarity,  the 
thermal  conductivity  of  the  mixture  is  larger  than  that  calculated 
by  an  averaging  procedure  (based  on  mole  fractions)  while  for  non-polar 
gas  mixtures  the  opposite  is  true.   The  difference  is  more,  the  more 
the  differences  in  the  sizes  of  the  different  constituent  molecules. 


189 


The  experimental  determination  of  the  thermal  conductivity 
of  mixtures  is  complicated  since  it  is  often  inseparably  tied  up 
with  other  phenomena.   Thus  the  presence  of  a  temperature  gradient 
causes  diffusion  (Soret  effect)  which  in  turn  sets  up  a  concentration 
gradient  resulting  in  the  transport  of  energy  (Dufour  effect);  the 
latter  is  generally  very  small  and  negligible,  even  when  the  former 
is  not. 

Many  authors  have  suggested  various  em.pirical  (for  particular 
mixtures),  semi -empirical  and  theoretical  formulae  for  the  determination 
of  the  mixture  thermal  conductivity:   Vasil'yeva,  Enskog,  Chapman, 
Lindsay  and  Bromley,  Mason  and  Saxena,  etc.;  an  extensive  compilation 
may  be  found  in  Refs.  18  and  27. 

For  the  special  case  of  a  binary  mixture,  the  assumption  of 
the  Lennard- Jones  parameters  for  the  entire  mixture,  as  those  obtained 
by  the  combining  rules  (2)  and  (3)  would  give  fairly  satisfactory 
results . 

Here  the  semi-em.pirical  formula  of  Wilke  [28]  has  been  used 
for  estimating  the  mixture  thermal  conductivity.   The  agreement  of 
this  formula  with  the  limited  experimental  data  is  fairly  good. 


N 


x.X. 
1  1 


mxx    .^^       N 


(21) 


where   d) .  .  = 


1  + 


M. 

1 

M. 


-1/2 


1  + 


,— r2 


^1/4 

(22). 
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This  is  the  only  equation  where  the  coefficient  of  viscosity 
appears  in  this  analysis;  hence  the  above  discussion  has  included 
the  determination  of  viscosity  also,  though  it  does  not  appear 
explicitly  in  the  conservation  equations. 

Transport  properties  of  free  radicals.   Experimental  data  on 
the  transport  properties  of  radicals  and  atoms  are  almost  non-existent. 
The  transport  phenomena  in  these  substances  are  masked  by  chemical 
reactions  and  related  phenomena,  posing  difficult  problems  to  get 
reliable  experimental  data.   Thus  one  has  to  again  rely  on  theoretical 
methods  to  get  any  estimate  at  all  of  the  transport  properties.   The 
theoretical  development  is  yet  in  a  rudimentary  stage  because  of  the 
complex  nature  of  the  phenomena.   Molecules  with  open  outer  electronic 
shells  exhibit  markedly  different  behavior  from  the  valence  saturated 
molecules.   The  Lennard-Jones  (6-12)  potential  cannot  be  used  to 
predict  the  transport  properties  of  atoms  and  radicals.   For  large 
distances  of  separation,  however,  the  behavior  is  similar  since  the 
electronic  wave  functions  of  the  two  molecules  essentially  do  not 
overlap,  so  that  there  is  the  usual  seventh  pov;er  van  der  Waals  force 
of  attraction.   As  the  tv/o  molecules  coma  together  there  are  several 
possible  interaction  potentials  corresponding  to  different  alignments 
of  the  electronic  spins,  some  attractive  and  the  others  repulsive. 
Using  quantum  mechanics  and  Pauli's  exclusion  principle  the 
probability  of  each  one  of  these  states  can  be  predicted  and  the 
problem  can  be  handled  by  suitably  weighted  averages  [29].   Yun  and 
Mason  [30]  have  utilized  this  m.ethod  to  construct  interaction  potentials 
for  the  binary  pairs  N-N   0-0   and  0-0,  and  the  transport  properties 
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for  these  pairs  based  on  a  set  of  collision  integrals  have  been 
calculated  [31,32],   In  the  absence  of  experimental  data,  verifica- 
tion of  this  theory  is  not  possible.   However,  there  is  fairly  good 
agreement  between  the  calculated  values  for  the  pair  0-0  extrapolated 
to  room  temperature  and  the  experimentally  measured  values. 

Some  simple  potentials  incorporating  an  inverse  po'-zer  repulsion 
law  of  the  form  k'T   (with  the  k'  and  s  determined  by  best  fits  to 
experimental  data,  if  available)  may  also  be  used  to  obtain  some 
estimates  of  the  transport  properties  [24]. 

In  some  flames,  some  of  the  species  may  be  in  electronically 
excited  states.   In  such  cases  it  may  be  best  to  consider  these  as 
totally  different  species  for  the  purpose  of  analysis.   Due  to  the 
large  collision  cross-section  of  such  species,  their  diffusion 
coefficients  would  essentially  be  small  [33]  and  may  be  neglected 
to  a  first  approximation.   The  transport  phenom.ena  would  also  be 
obscured  by  electronic  exchanges  and  chemical  reactions,  making 
theoretical  and  experimental  determination  of  the  properties  at 
best  difficult. 

In  the  analysis  carried  out  here,  three  different  types  of 
atoms  (H,  0,  N) ,  and  three  radicals  (CH^ ,  CHO,  OH)  have  been 
considered  and  they  have  been  assumed  to  be  in  the  electronic 
ground  state,  which  is  a  good  assumption  at  the  temperature  existing 
in  this  flame  (approximately  2400 °K).   The  theoretical  predictions 
of  Refs.  31  and  32  for  H-H   N-N   and  0-0   (reproduced  below)  have 
been  used  to  fit  to  a  Lennard-Jones  (6-12)  potential  and  calculate 
the  parameters  O   and  £  for  these  species.   This  was  done  in  the 
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absence  of  other  suitable  techniques  to  evaluate  the  transport 
properties,  realizing  fully  well  that  this  type  of  interaction 
potential  was  not  valid  here.   Moreover  the  parameters  O   and  £ 
cease  to  have  their  usual  physical  meaning  in  this  context;  however, 
they  may  be  regarded  as  hypothetical  quantities  which  when  used  in 
conjunction  with  the  Lennard-Jones  (5-12)  potential  would  give 
estimates  of  the  transport  properties  which  agree  reasonably  well 
with  the  limited  experimental  data  available.   For  the  species  CH 
CHO,  and  OH,"  since  no  data  whatsoever  could  be  found  in  the 
literature,  'suitable'  values  had  to  be  guessed  and  used. 

Since  all  these  species  are  present  in  very  small  concentra- 
tions relative  to  the  major  species,  even  large  errors  in  the 
estimation  of  their  transport  properties  would  affect  the  overall 
analysis  only  slightly. 

Table  1 
Atom-molecule  diffusion  coefficients  computed  from  theory 

^Ml_Y§-'-!i£§  at  1  atm.  pressure) 

2 
Diffusion  coefficient  (en,  /sec.) 


Temperature 

H-H^ 

N-N^ 

0-0^ 

(degrees 

K) 

Ref.    31 

Ref.    32 

Ref.    32 

1000 

17.2 

2.34 

2.40 

1500 

36.3 

4.78 

4.80 

2000 

62.0 

7.95 

7.87 

2500 

94.2 

11.8 

11.6 

3000 

133.0 

16.4 

15.9 

3500 

178.0 

21.6 

20.7 

4000 

230.0 

27.5 

26.2 

■k 
For  OH,  Wilde  [34]  has  used  a  certain  value  of  a  dimensionless 
diffusion  parameter;  the  source  of  the  data  is,  however,  not  quoted. 
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Thermal  diffusion.   Thermal  diffusion  is  not  of  much 
significance  in  most  flames  and  is  consequently  neglected.   Its 
contribution,  if  any,  would  be  towards  the  diffusion  of  light  gases 
such  as  atomic  and  molecular  hydrogen.   Experimental  measurements 
of  the  thermal  diffusion  ratio  for  some  gas  pairs  at  low  temper- 
atures have  been  reported  [24,  p.  584].   To  obtain  values  of  this 
quantity  at  flame  temperatures,  theoretical  calculations  are  the 
only  available  alternative.   However,  this  quantity  is  m^ore  sensitive 
to  the  form  of  interaction  potential  chosen  than  the  other  transport 
properties  and  so  theoretical  calculations,  especially  at  temperatures 
encountered  in  flames,  are  at  best  approximate. 

If  it  is  desired  to  consider  the  thermal  diffusion  for  any 
or  all  of  the  species,  it  is  convenient  to  treat  the  given  multi- 
component  mixture  as  a  binary  mixture,  one  component  of  which  is  the 
particular  species  being  studied,  the  other  species  being  all  lumped 
into  the  second  component.   Then  the  contribution  to  the  diffusion 
velocity  due  to  thermal  diffusion  is  proportional  to  (k  /T)* (dT/dz) , 
where  k   is  the  thermal  diffusion  ratio  for  the  hypothetical  binary 
mixture.   An  order  of  magnitude  analysis  reveals  that  this  contribu- 
tion to  diffusion  velocity  is  generally  an  order  of  magnitude  lov/ar 
than  the  contribution  by  the  concentration  gradient,  so  its  neglect 
is  justified  [35,  p.  50], 

Dif fusion-thermo,  the  transfer  of  energy  due  to  a  concentration 
gradient  is  generally  negligible  even  v;hen  thermal  diffusion  is  not 
negligib]e . 
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Limitations  of  the  Chapman-Enskog  theory.   The  determination 
of  the  transport  properties  of  viscosity,  thermal  conductivity  and 
diffusion  as  outlined  above  is  largely  dependent  on  the  Chapman-Enskog 
theory.   Therefore,  before  concluding,  it  would  not  be  out  of  place 
to  mention  some  of  the  limitations  of  this  theory. 

(i)   The  Chapman-Enskog  method  for  the  solution  of  the 
Boltzmann  equation  yields  a  series  solution,  only  the  first  term 
of  which  has  been  considered  here.   The  first  order  theory  is  valid 
when  distances  over  which  appreciable  changes  in  macroscopic  properties 
(like  pressure,  temperature,  density,  etc.)  occur  are  large  compared 
to  a  mean  free  path.   In  mathematical  terms  this  means  £  (8ilnl|^/8x)  «  1 
where  i   =  m£an  free  path,  and  ^  =   any  macroscopic  property.   Since 
the  m.ean  free  path  is  inversely  proportional  to  the  pressure,  for  a 
given  large  gradient  in  a  physical  property,  the  theory  would  break 
down  sooner  at  lower  pressures.   The  first  order  theory  considers 
the  fluxes  being  proportional  to  the  first  derivative  of  the 
corresponding  physical  properties.   The  higher  order  terms  involve 
higher  derivatives  and  non-linear  terms  involving  powers  and  products 
of  lower  derivatives. 

(ii)   The  theory  considers  only  binary  collisions,  hence  the 
densities  should  be  low  enough  that  this  approximation  be  valid. 
Tlie  equations  for  the  transport  coefficients  based  on  the  theory  are 
generally  valid  for  pressures  less  than  about  20  atmospheres.   The 
equations  fail  to  give  correct  results  at  high  pressure,  since, 
besides  considering  only  binary  collisions,  the  theory  does  not 
account  for  the  effect  of  molecular  size  on  the  frequency  of 
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molecular  collision,  which  is  important  at  high  pressures.   The 
theory  of  Enskog  gives  fairly  reasonable  estimates  of  the  transport 
properties  of  dense  gases.   Other  more  rigorous  theories  are  avail- 
able, but  they  are  not  as  convenient  to  use.   Hirschfelder  et  al. 
[24]  have  given  some  results  to  indicate  the  agreement  of  Enskog's 
theory  with  experimental  data  for  the  coefficient  of  viscosity. 
The  development  of  this  theory  for  the  diffusion  coefficient  has 
been  less  successful.   There  is  also  a  lack  of  experimental  data 
at  high  pressures  so  that  the  accuracy  of  the  theoretical  predictions 
cannot  be  checked. 

(iii)   The  theory  assumes  that  molecular  interaction  is  described 
by  classical  mechanics.   Quantum  effects  might  become  prominent  at  low 
temperatures  and  high  densities,  but  are  almost  always  negligible  in 
ordinary  flame  studies. 

(iv)   The  dimensions  of  the  container  and  any  obstacles  within 
it  are  assumed  to  be  large  enough  compared  to  the  mean  free  path,  so 
that  the  volume  of  the  surface  layer  is  negligible  compared  to  the 
total  system  volume. 

(v)   The  theory  applies  strictly  only  to  monatomic  gases;  i.e., 
molecules  which  have  no  internal  degree  of  freedom  and  have  a  spherically 
symmetric  interaction  potential. 

(vi)   Besides  the  physical  limitations  m.entioned  above,  the 
Chapman-Enskog  theory  has  a  serious  drawback  in  its  application. 
Once  an  interaction  potential  model  is  chosen,  the  calculations 
required  are  considerable,  so  that  the  theory  permits  only  the 
determination  of  the  best  parameters  by  comparison  with  experimental 
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data  in  a  previously  chosen  model,  rather  than  the  determination  of 
the  best  interaction  model.   This  causes  the  parameters  of  a  given 
potential  to  be  different  when  evaluated  by  best  fits  to  experimental 
data  of  different  properties ;  again  two  different  models  giving  the 
same  final  result  may  yield  significantly  different  interaction 
potentials . 
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APPENDIX  III 
COMPUTER  PROGRAIIS 

In  this  appendix,  a  listing  of  the  computer  programs  used 
in  this  work  is  provided.   Sufficient  explanation  is  included  in 
the  programs  so  that,  together  with  the  theory  developed  in  the 
preceding  chapters,  they  are  more  oi"  less  self-explanatory. 
Tlie  programs  have  been  split  into  three  parts:   the  first  one  deals 
with  the  analysis  of  the  Otto  cycle  from  the  point  of  fuel-air  irdxture 
intake  to  the  point  where  combustion  is  completed;  the  second  one 
deals  with  the  expansion  stroke  analysis  using  the  results  from  the 
first  program;  and  finally  the  third  program  deals  with  the  structure 
of  the  flame  front  and  its  propagation. 

For  each  one  of  the  programs  the  parameters  of  interest, 
memory  space  required,  and  execution  time  are  specified.   The 
prograffiS  are  all  written  in  FORTRAN  IV  and  the  execution  times 
indicated  are  those  obtained  when  using  the  FORTRAN  H  EXTENTDED 
(0PT=2)  compiler  on  an  IBM  370/165  computer.   Some  representative 
output  is  also  included. 

In  all  three  programs  the  thermodynamic  data  have  been 
obtained  from  the  JANAF  Tables  by  storing  the  values  of  the  prop- 
erties of  interest  over  the  temperature  range  likely  to  be  encountered. 
This  requires  considerable  memory  storage  space  and  input  data.   An 
alternate  but  less  precise  method  would  be  to  utilize  polynomial  fits 
by  the  method  described  in  Ref.  1.   A  similar  remark  applies  to  the 
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table  of  collision  integrals  used  in  the  evaluation  of  transport 
properties  in  Program  III. 


P   R   0  C-   R  A  M I 

This   program  is    used   for   the   analysis   of   the   Otto   cycle 
from  the   time   of   intake   of   the    fuel-air  mixture    until   the    time 
combustion   is    completed.       Subroutine    SE-IQ   used    for   the    simultaneous 
solution   of   a  set   of   algebraic  equations   has   been   taken    from  the 
IBM   Scientific    Subroutine   Package. 

Memory   space   required:       75,000  bytes 

Execution   time:      9   seconds 
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PROGRAM      II 

This   program  is    used   for   the   analysis    of   the   expansion 
stroke    from  the    time    combustion   is    completed   until   exhaust. 
Results    from  program   I  are    used   as    input    data.       Subroutines 
DECOMP,    IMPRUV,    and   SOLVE  have  been   taken    from  Ref.    2.       Subroutine 
REACT   from  Ref.    3  has  been  modified   extensively  while    retaining 
the   same  basic   structure   so   as    to  be    of  more   relevance   to   the 
problem  under  study  here. 

Memory   space   required:      92,000  bytes 

Execution   time:      3  minutes 
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This  program  is  used  for  the  numerical  solution  of  the  set 
of  partial  differential  equations  governing  one-dimensional  unsteady 
laminar  prem.ixed  flame  propagation.   Subroutines  DECOMP,  F,  F2 , 
E-IPRUV,  PEDERV,  REACT,  SOLVE  are  identical  with  those  used  in 
Program  II  and  are  not  listed  here.   All  these  subroutines  are 
involved  in  the  integration  of  a  set  of  stiff  ordinar:>'  differential 
equations.   Subroutines  SIMQ  and  QSF  are  from  the  IBM  Scientific 
Subroutine  Package. 

Memory  space  required:   222,000  bytes 

Execution  time:   This  depends  strongly  on  the  starting  profiles  used 
for  temperature  and  concentrations.   With  the  profiles  used  here,  the 
time  is  approximately  20  minutes. 
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